Introducción

El análisis de series climáticas es siempre un desafío y más aún en regiones donde la disponibilidad de las mismas es escasa. Este déficit es un problema importante para la realización de análisis de riesgo o de modelación que demandan una gran cantidad de datos. Los generadores estocásticos de clima producen series diarias de variables climáticas con propiedades estadísticas consistentes a los de los datos observados. Estas series sintéticas son un insumo muy importante para el análisis cuantitativo del riesgo de desastres naturales. Este tipo de aplicaciones requiere de largas series temporales para poder capturar la variabilidad natural del clima de una región. Las principales variables de importancia son temperaturas máxima y mínima y precipitación. En muchas regiones estos datos se encuentran incompletos o son inexistentes. Los registros suelen tener una duración insuficiente o sólo estar disponibles datos agregados mensualmente. Para suplir este déficit, los generadores estocásticos pueden generar series que satisfagan las necesidades.

La mayoría de los enfoques tradicionales para la generación de series estocásticas están limitados por su capacidad de generar datos solamente en localidades para las que se cuenta con observaciones (por ejemplo, donde existen estaciones meteorológicas). Otra desventaja de algunas de estas herramientas (sobre todo los generadores no paramétricos basados en remuestreo de observaciones) es que solamente pueden producir valores dentro del rango observado en el registro histórico. En este proyecto, el generador desarrollado se basó en el modelo diseñado por Verdin et al. (2016). Este generador estocástico diario multivariado produce series sintéticas de precipitación y temperaturas máxima y mínima. El generador de Verdin et al. (2016) utiliza cuatro modelos estadísticos para modelar la ocurrencia y montos de precipitación y temperatura máxima y mínima basados en Modelos Lineales Generalizados (GLMs). Estos modelos son paramétricos y utilizan regresiones lineales entre las variables para modelarlas. El proceso de ocurrencia de precipitación se modela como una regresión probit mientras que los montos se ajustan a una distribución aleatoria Gamma. Las temperaturas máximas y mínimas son modeladas como variables autorregresivas condicionadas por la precipitación. El generador produce series espacialmente correlacionadas tanto para estaciones individuales como en grillas regulares. Una de las limitaciones fundamentales del generador de Verdin et al. (2016) es que las temperaturas máxima y mínima se generan en forma independiente para cada día, por lo cual la amplitud térmica diaria simulada a veces no es realista. Los diagnósticos previos realizados en base al generador de Verdin et al. (2016) demostraron que algunas propiedades de las series sintéticas producidas no reflejaban fielmente características importantes para el análisis de las sequías, por ejemplo, la persistencia de secuencias de días secos y lluviosos. Por este motivo, en este trabajo solo se mantuvo la estructura general del modelo de Verdin et al. (2016) y se modificaron todos los algoritmos (modelos estadísticos) para obtener series sintéticas más consistentes con los registros históricos y para optimizar el proceso de cálculo.

Fundamento

A partir de los datos climáticos históricos, comienza el proceso de ajuste de un modelo estadístico que permita representar el comportamiento de cada una de las variables para cada localidad. El generador está dividido en cuatro modelos aditivos generalizados: dos para modelar la precipitación y dos para las temperaturas máxima y mínima, respectivamente. El concepto principal detrás de la generación de series sintéticas es que cada valor de una variable puede ser considerado como la suma de una componente climática (clima local) y otra componente meteorológica aleatoria o tiempo local (Kleiber et al., 2013), es decir

\[ X_{i,s} = clima\ local + tiempo\ local \] donde \(X_{i,s}\) corresponde al valor de la variable X en el día i en la localidad s; clima local corresponde a un valor medio de la variable para el día i en la localidad s y tiempo local corresponde a un estado particular de la atmósfera en el día i en la localidad s. El componente climático es especificado a través del ajuste de cada uno de los cuatro modelos aditivos del generador. El componente meteorológico corresponde a los residuos de los modelos (la diferencia entre los valores históricos y el componente climático estimado), es decir, a la variabilidad no explicada por los mismos. El proceso de ajuste culmina con los parámetros ajustados para los cuatro modelos mencionados. Posteriormente, se generan datos para los años a simular que corresponden a los valores medios de la distribución para cada día del año (clima local). Luego, se simulan una serie de valores aleatorios, a partir de los residuos de cada modelo, que corresponden a la variabilidad propia de cada realización (tiempo local).

Figura 1.  Componentes clima local y tiempo local

Figura 1. Componentes clima local y tiempo local

Tipos de series sintéticas

El generador GAMWGEN produce distintos tipos de series sintéticas. Desde el punto de vista temporal, las series pueden ser no condicionadas (Figura 2.a), es decir, puramente estacionarias; pseudohistóricas (Figura 2.b), que copian la variabilidad de baja frecuencia y los cambios en la serie climática observada y condicionadas (Figura 2.c) que son series forzadas a seguir la trayectoria de una salida de un modelo de cambio climático o una trayectoria arbitraria definida por el usuario. Las siguientes Figuras muestran distintos ejemplos de lo mencionado.

Figura 2.a.  Tipos de series sintéticas: Series estacionarias

Figura 2.a. Tipos de series sintéticas: Series estacionarias

La línea negra corresponde a la serie temporal observada mientras que las distintas líneas celestes corresponden a realizaciones del generador. Al tratarse de series puramente estocásticas, son independientes entre sí.

Figura 2.b.  Tipos de series sintéticas: Series pseudihistóricas

Figura 2.b. Tipos de series sintéticas: Series pseudihistóricas

La línea negra corresponde a la serie temporal observada mientras que las distintas líneas celestes corresponden a realizaciones del generador. En este caso, el generador fue forzado a seguir la trayectoria en los datos observados, por lo tanto, las series generadas siguen las variaciones de los datos observados.

Figura 2.c.  Tipos de series sintéticas: Series condicionadas

Figura 2.c. Tipos de series sintéticas: Series condicionadas

En esta configuración el generador sigue una trayectoria arbitraria elegida por el usuario. En este caso es una trayectoria lineal por lo que la precipitación aumenta un determinado porcentaje por año de manera lineal.

Desde el punto de vista espacial se pueden generar series en estaciones meteorológicas, en puntos que no se correspondan con estaciones o en una grilla regular. La Figura 3 muestra un ejemplo de la generación en grillas sobre el Paraguay.

Figura 3.  Tipos de series sintéticas: Datos grillados

Figura 3. Tipos de series sintéticas: Datos grillados

El panel de la izquierda muestra la temperatura máxima del día 1 de enero de 2000 para todo el territorio del Paraguay mientras que el mapa del panel derecho muestra la precipitación diaria para el mismo día.

Metodología

Como su nombre indica, el generador está basado en Modelos Generalizados Aditivos (GAM, por sus siglas en inglés). Como se mencionó en le Introducción, este generador está basado en uno similar desarrollado por Verdin et al. (2016). Dicho generador estocástico utilizaba modelos lineales generalizados (GLM, por sus siglas en inglés). Los GLM son modelos muy interesantes ya que son facilmente interpretables aunque carecen de la flexibilidad necesaria para capturar complejos patrones como las variaciones estacionales. Por ello, en esta versión se cambiaron los GLMs por GAMS. Estos nuevos modelos están siendo muy utilizados en diversas áreas porque heredan lo mejor de los GLM pero son mucho más flexibles. En la siguiente sección se describirán brevemente los GAMs.

Introducción a los Modelos Aditivos Generalizados

Interpretabilidad vs Complejidad

Es posible crear modelos ajustados a los datos que son lineales y relativamente fáciles de interpretar y explicar. Existe una línea recta que representa la relación entre dos variables y atraviesa la nube de puntos. Es posible graficar la relación y hacer predicciones a partir del modelo ajustado. Pero lo modelos lineales no siempre representan bien las relaciones entre variables ya que las relaciones no siempre son lineales. Las predicciones no serían buenas a partir de estos modelos.

En el otro extremo del espectro hay toda una serie de modelos de tipo “caja negra” como las redes neuronales, random forest, árboles de regresión. Estos modelos son muy buenos para predicir pero son muy díficiles de interpretar y de entender que está sucediendo en el sistema. Son muy útiles para clasificar pero no sirven para entender como una variable se relaciona con el producto del modelo.

Los GAMs proveen un interesante punto intermedio ya que se pueden ajustar relaciones complejas y no lineales e interacciones pero estos modelos son explicitos y se pueden observar las relaciones entre variables y entender porque se produce un resultado determinado.

Relaciones no lineales

En general, las relaciones entre variables de la naturaleza no son lineales y adquieren patrones muy complejos. En la Figura 4 se muestra un ejemplo de datos sintéticas para ejemplificar este concepto.

library(mgcv)
set.seed(2) ## simulate some data...
dat <- gamSim(1, n=400, dist="normal", scale=0, verbose=FALSE)
dat <- dat[,c("y", "x0", "x1", "x2", "x3")]
ggplot2::ggplot(dat, ggplot2::aes(y=y,x=x2)) +
    ggplot2::geom_point() +
    ggplot2::theme_minimal()
Figura 4.  GAMs: Relaciones no lineales.

Figura 4. GAMs: Relaciones no lineales.

Este diagrama de dispersión muestra dos variables que claramente se encuentran relacionadas pero no de manera lineal. Si se ajusta una regresión lineal simple a estos datos se obtienen los siguientes resultados.

ggplot2::ggplot(dat, ggplot2::aes(x2, y)) + 
    ggplot2::geom_point() +
    ggplot2::geom_smooth(method = 'lm', se = TRUE, ggplot2::aes(colour = "Lineal")) +
    ggplot2::scale_colour_manual(name = "", values = c("Steelblue")) +
    ggplot2::theme_minimal()
Figura 5.  GAMs: Relaciones no lineales. Ajuste lineal. Tomado de (CITAR)

Figura 5. GAMs: Relaciones no lineales. Ajuste lineal. Tomado de (CITAR)

El modelo no es capaz de capturar las principales características de los datos.

Modelo lineal
modelo_lineal <- lm(y ~ x2, data = dat)

Gráfico de residuos

residuals <- fortify(modelo_lineal)
ggplot2::ggplot(residuals, aes(x = .fitted, y = .resid)) + 
  ggplot2::geom_point() +
  ggplot2::geom_smooth(se = FALSE) +
  ggplot2::geom_hline(yintercept = 0, linetype = 'dashed') +
  ggplot2::theme_bw() +
  ggplot2::xlab('Ajustados') + ggplot2::ylab('Residuos')
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Claramente los residuos tienen una tendencia y no están distribuidos aleatoriamente.

Sin embargo, al cambiar el modelo lineal por un gam.

Figura 6.  GAMs: Relaciones no lineales. Ajuste no lineal. Tomado de (CITAR)

Figura 6. GAMs: Relaciones no lineales. Ajuste no lineal. Tomado de (CITAR)

Con un GAM se pueden ajustar los modelos a través de funciones suavizadas o splines que pueden tomar casi cualquier forma. Al utilizar splines los GAMs pueden capturar diversos tipos de relaciones no lineales y es por esto que son tan flexibles y adaptables a distintos contextos.

modelo_gam <- gam(y ~ s(x2), data = dat)
residuals <- data.frame(fitted = modelo_gam$fitted.values,
                            resid = resid(modelo_gam)) %>%
    tibble::as_tibble()
ggplot2::ggplot(residuals, aes(x = fitted, y = resid)) + 
  ggplot2::geom_point() +
  ggplot2::geom_smooth(se = FALSE) +
  ggplot2::geom_hline(yintercept = 0, linetype = 'dashed') +
  ggplot2::theme_bw() +
  ggplot2::xlab('Ajustados') + ggplot2::ylab('Residuos')
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Al utilizar un GAM los residuos tienen una distribución casi aleatoria sin una tendencia clara como si fue el caso del modelo lineal. Este es sólo en sencillop ejemplo del potencial que tienen los GAMs para modelar complejas relaciones. Si se desea profundizar más en el tema se recomienda revisar el librop de Simon Wood (2017)

Modelación

Como se mencionó anteriormente, el generador estocástico tiene como base cuatro modelos generalizados, dos para temperaturas máxima y mínima y dos para la precipitación que modelan la ocurrencia y los montons diarios.

Figura 6.  Modelos estadísticos

Figura 6. Modelos estadísticos

Partiendo de la base de datos se extraen los datos diarios de temperaturas máxima y mínima y precipitación. Los valores de temperaturas son utilizados para ajustar un GAM que modelará la componente climática. Estos modelo puede incluir o no medias trimestrales que servirán para condicionar el modelo a seguir un determinado comportamiento. Para el caso de la precipitación el tratamiento es diferente. El fenómeno se divide en dos: ocurrencia de precipitación y montos diarios. La ocurrencia se modela con GAM para toda la serie que puede o no estar condicionado por totales trimestrales de lluvia. Los montos, en cambio, son modelados a escala mensual. Es decir, se ajustan doce modelos, uno para cada mes del año. Se utiliza esta modalidad para capturar mejor las características propias del ciclo estacional de la precipitación de cada región.

A continuación se describirán cada uno de los modelos recién mencionados.

Modelación: Clima local: Ocurrencia de lluvia

El ajuste del modelo de ocurrencia comienza con la definición de un día lluvioso. Según la Organización Meteorológica Mundial (OMM), se considera como día lluvioso a aquel día con una precipitación igual o mayor a 0.1 mm. Este valor, si bien es el sugerido por la OMM, puede ser modificado por el usuario si así lo dispone. Una vez definida la ocurrencia de lluvia se ajusta el modelo. La ocurrencia de lluvia es una variable de tipo binaria, es decir, un día puede ser lluvioso (1) o seco (0), y es muy bien representada a través de una regresión probit Kleiber et al. (2012). Este tipo de regresiones se basan en procesos latentes Gaussianos, \(W_{s,t}\), que se modelan con la siguiente relación:

\[ O_{s,t} = \Pi_{\{W_{s,t} > 0|} \] Si el proceso \(W_{s,t}\) es positivo, significa que lloverá en el día t y en la estación s y a ese día se le asignará el valor 1. Si el proceso es negativo significa que no lloverá en el día t y en la estación s y ese día se le asignará el valor 0. El uso de este tipo de procesos latentes está justificado en que, en una región, la ocurrencia de lluvia en las distintas estaciones tenderá a estar correlacionada. La función media del proceso latente Gaussiano es una regresión entre variables que se expresa de la siguiente manera:

\[ O_{s,t} = (s, O_{s,t-1}, f(doy(t)), f(ST(t)), f(lon, lat)) \] donde \(O_{s,t}\) corresponde a la ocurrencia de lluvia en sitio y día determinado; \(s\) corresponde al efecto del sitio s (ordenada al origen); \(O_{s,t-1}\) corresponde a la ocurrencia del día previo que es un término autorregresivo; \(f(doy(t))\) corresponde a una función cíclica de los días del año para considerar el efecto de la estacionalidad sobre la ocurrencia de lluvia; \(f(ST(t))\) corresponde a una función suavizada de los acumulados estacionales de precipitación (solo se utiliza si se desea condicionar el modelo) y \(f(lon, lat)\) corresponde a las coordenadas del punto s (solo se utiliza en el modelo espacial). En la práctica, esta covariable se divide en cuatro, una para cada trimestre del año, asignándole el valor de 0 para los momentos fuera del respectivo trimestre. Es importante notar que para la ocurrencia de precipitación se usa un solo término autorregresivo. Este término es muy importante para la correcta modelización de las rachas secas y lluviosas mientras que el día del año incorpora la variabilidad intra-anual.

Modelación: Clima local: Montos diarios de lluvia

El modelo de montos diarios de lluvia difiere del anterior en que no se usan todos los datos de la serie, sino que se extraen de la base de datos los montos de precipitación únicamente para los días lluviosos. La intensidad de la precipitación para una localidad s y tiempo t es modelada como una variable aleatoria Gamma cuyos parámetros de forma y escala varían en el tiempo y en el espacio (Kleiber et al. (2012)). La función Gamma ha sido ampliamente utilizada en la región para la modelación de acumulados de lluvia. El modelo de montos puede ser expresado de la siguiente manera:

\[ I_{s,t} = (s, O_{s,t-1}, f(ST(t)), f(lon, lat)) \]

donde, \(I_{s,t}\) corresponde a los montos de precipitación en un sitio y día determinado; \(O_{s,t-1}\) corresponde a la ocurrencia del día previo para considerar la autocorrelación temporal; \(f(ST(t))\) corresponde a una función suavizada de los acumulados estacionales de precipitación (solo se utiliza si se desea condicionar el modelo) y \(f(lon, lat)\) corresponde a las coordenadas del punto s (solo se utiliza en el modelo espacial). A diferencia del modelo anterior – que modela los años calendarios completos a través de la estacionalidad del día del año – la serie de montos diarios es dividida en función del mes del año para ajustar una distribución a cada mes para obtener así parámetros mensuales más precisos. Gracias a esta modificación se obtienen parámetros que varían en el tiempo y en el espacio lo que permite capturar la variabilidad espacial de la precipitación. Al igual que en el modelo de ocurrencia, la inclusión del total trimestral es necesaria si se desean simular series condicionadas.

Modelación: Clima local: Temperatura

Para el caso de la temperatura se utiliza una metodología similar a la utilizada para la precipitación, basada en Kleiber et al. (2013). A partir de los datos observados de temperatura se ajustan dos modelos: uno de máxima y otro de mínima. Ambos modelos utilizan las mismas variables para realizar el ajuste. El modelo puede ser expresado de la siguiente manera:

\[ X_{s,t} = (s, O_{s,t}, O_{s,t-1}, f(T_{x_{(s, t-1)}}, T_{n_{(s, t-1)}}), f(doy(t)), f(SX(t), SN(t)), f(lon, lat)) \] donde, \(X_{s,t}\) corresponde a la temperatura máxima o mínima en un sitio y momento determinado; \(s\) corresponde al efecto del sitio s (ordenada al origen); \(O_{s,t}\) corresponde a la ocurrencia de lluvia en sitio y día determinado; \(O_{s,t-1}\) corresponde a la ocurrencia del día previo; \(f(T_{x_{(s, t-1)}}, T_{n_{(s, t-1)}})\) corresponde a una función suavizada de la interacción entre la temperatura máxima y mínima del día previo; \(f(doy(t))\) corresponde a una función cíclica de los días del año para considerar el efecto de la estacionalidad sobre la temperatura; El término \(f(SX(t), SN(t))\) corresponde a la interacción entre las medias estacionales de temperatura máxima y mínima y, como se explicó en la sección anterior, se incluyen para crear modelos condicionados y \(f(lon, lat)\) corresponde a las coordenadas del punto s (solo se utiliza en el modelo espacial). La ocurrencia de lluvia es muy importante ya que en general los días lluviosos tienen temperaturas más bajas que los secos, sobre todo en verano, por lo que debe ser incluido en el ajuste. La Figura 6 es un ejemplo de esta influencia en Junín, en donde se muestra la amplitud térmica para cada mes del año en función del tipo de día, seco o lluvioso.

Modelación: Tiempo local: Para una estación

Precipitación

El tiempo local de la ocurrencia de precipitación se modela a través de los residuos de la regresión probit. Por definición estos residuos tienen una distribución \(X \sim \mathcal{N}(0,\,1)\,\). Por lo tanto, generar valores para el tiempo local es muy sencillo, solo se necesita de una función gaussiana que genere números aleatorios. La Figura muestra un ejemplo de la generación de clima local y tiempo local para una estación meteorológica de Argentina.

Figura 6.  Tiempo local: precipitación

Figura 6. Tiempo local: precipitación

Esta Figura muestra sólo un año para mejorar la visualización pero su interpretación es válida para toda la longitud de la serie. La línea azul corresponde al clima local modelado a través del GAM y la naranja al ruido aleatorio creado a partir de una distribución \(X \sim \mathcal{N}(0,\,1)\,\). La sumatoria de ambos componentes determinarán si el día es lluvioso o no. Si la suma es positiva, lloverá, caso contrario será un día seco. Se puede observar en la línea azul un patrón estacional debido al régimen de precipitación tipo monzónico de esta región con picos de más días lluviosos durante el verano mientras que en invierno disminuyen marcadamente. Esta misma serie temporal de números aleatorios serán la base para la generación de los montos de precipitación.

Temperatura

El tiempo local de las temperaturas máxima y mínima se modela de una manera diferente. En este caso se toman los residuos de cada uno de los GAMs, es decir, la diferencia entre el valor ajustado por el modelo y el valor observado de temperatura. Además, como la temperatura está fuertemente influenciada por el tipo de día, días lluviosos tienden a tener una menor amplitud térmica que los días secos, los residuos se separan en función del tipo de día. Para capturar mejor el patrón estacional de la temperatura, los residuos se agrupan por mes y se ajusta un modelo bivariado que contemple los residuos de temperaturas máxima y mínima. El uso de un modelos bivariado es una alternativa muy interesante para mantener la consistencia entre ambas variables, es decir, que la temperatura máxima no supere a la mínima. La siguiente Figura es un ejemplo de las series de tiempo local para una localidad de Argentina.

Figura 7.  Tiempo local: temperatura

Figura 7. Tiempo local: temperatura

La figura esta dividida en dos paneles, el superior muestra el tiempo local para los días lluviosos y el inferior para los días secos. La línea naranja corresponde a los valores de tiempo local para temperatura máxima y los azules para temperatura mínima. Si bien ambas series difieren en magnitud y variabilidad, las dos tienden a variar conjuntamente.

Modelación: Tiempo local: Para una grilla

Para generar datos sobre una grilla regular o en puntos donde no hay datos para ajustar los modelos se debe utilizar el modelo espacial. La generación del tiempo local en el espacio es conceptualmente idéntico a la generación sobre estaciones meteorológicas sólo que utiliza campos gaussianos para incluir la dimensión espcial. Los campos gaussianos se generan con el paquete RandomFields y capturan la variabilidad espacial de cada una de las variables a partir de su variograma.

Precipitación

Para la precipitación se utiliza un modelo exponencial que utiliza como parámetros el variograma ajustado a partir de los datos observados usando máxima verosimilitud. Este modelo permite simular campos con una muy buena consistencia espacial. La ocurrencia de precipitación no ocurre de manera aislada en una región sino que, en general, un evento lluvioso abarca una importante superficie. La siguiente Figura es una ejemplo de los campos aleatorios generados sobre una grilla regular para una región de Argentina.

Figura 8.  Tiempo local: precipitación sobre una grilla

Figura 8. Tiempo local: precipitación sobre una grilla

La interpretación es análoga a la mostrada para una estación puntual. Los valores para cada píxel se suman a la componente climática y así se determina si el día será lluvioso o no. Este tipo de campos se generan para todos los días de la simulación y cada campo diario es independiente del anterior.

Temperatura

Al igual que para una estación, los campos gaussianos que se generan son bivariados. Se utiliza el modelo Bivariado de Whittle Matern incluido en el paquete RandomFields. La siguiente Figura es un ejemplo para un día de una realización para grilla regular sobre Argentina.

Figura 8.  Tiempo local: temperatura sobre una grilla

Figura 8. Tiempo local: temperatura sobre una grilla

Ambos campos se generan simultáneamente para temperaturas máxima y minima y discriminando entre días secos y lluviosos. Esto quiere decir que en el proceso de creación de los campos se generan cuatro capas diferentes que luego se combinan en función del tipo de día de cada píxel resultando en dos campos integradores. Los valores para cada pi’xel se sumen a la componente climática y así se obtiene el valor final para cada día de la simulación.

Aplicación

En esta sección se mostrarán ejemplos de aplicación del generador para generar distintos tipos de series explicando las funciones necesarias y cada uno de los parámetros.

Instalar paquetes necesarios

El primer paso es comprobar que todos los paquetes necesarios estén instalados y si no es así, descargarlos e intalarlos.

## Loading required package: gamwgen

Creación de directorios

El paquete tiene precargados algunos ejemplos de aplicación con datos reales de estaciones meteorológicas de la red del SISSA.

Para poder seguir este tutorial se deben crear directorios donde se guardarán los datos de entrada y salida.

  • /input_data: aquí se guardarán los datos meteorológicos y los metadatos de las estaciones
  • /output_data: aquí se guardarán los resultados de la simulación

Si estos directorios no existen, se crearán.

Generación de series sintéticas para una sola estación meteorológica

Este primer ejemplo consiste en la generación de series sintéticas para estaciones meteorológicas. Es decir, se generarán series para las mismas estaciones que fueron utilizadas en el ajuste de los distintos modelos. Por lo tanto, no se considera la componente espacial. El ejemplo está divido en dos, en una primera parte se ajustará el modelo para generar series estacionarias y luego, en una segunda parte, se incluirán coviariables estacionales para producir series pseudohistóricas. Para ambos ejemplos los datos son usados serán los mismos.

Crear archivos de entrada

El primer paso consiste en generar los set de datos de entrada que se descargaron al momento de instalar el paquete del generador estocástico. Estos datos son sólo a título demostrativo, si el usuario desea correr el modelo con sus propios datos deberá cambiar los objetos que se generarán en esta sección por los suyos y colocarlos en la carpeta input_data.

Los archivos necesarios son:

  • stations.csv
  • climate.csv

Los datos meteorológicos se dividen en dos archivos separados: stations.csv y climate.csv. Los nombres de los mismos no deben ser necesariamente iguales a los usados aquí.

Los metadatos de las estaciones se alojan en el archivo stations.csv. Este archivo contiene la información de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:

  • station_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integer
  • latitude: latitud en grados decimales. La variable debe ser de tipo double
  • longitude: longitud en grados decimales. La variable debe ser de tipo double

La tabla puede tener más variables pero sólo se necesitan las anteriores.

A continuación se muestran la primera fila del dataset y los tipo de datos de cada una de las variables.

# Vista de los metadatos de la estación
knitr::kable(stations[1,])
x y station_id nombre lat_dec lon_dec elev pais_id
5001614 6256841 87448 Villa Reynolds Aero -33.7181 -65.3737 486 AR

El objeto stations debe ser convertido de tibble a sf. El sistema de referencia espacial debe ser planar. No es necesario un sistema de referencia espacial en particular, solamente las coordenadas deben estar expresadas en metros.

# Convertimos el objeto stations a sf y se convierte su proyección de WGS 1984 a
# POSGAR Argentina Faja 5.
stations %<>% 
  sf::st_as_sf(coords = c("lon_dec", "lat_dec"), crs =  4326) %>%
  sf::st_transform(crs = 22185)

La información climática se aloja en el archivo climate.csv. Este archivo contiene los datos de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:

  • date: fecha del dato. La variable debe ser de tipo date
  • station_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integer
  • prcp: datos diarios de precipitación La variable debe ser de tipo double
  • tmax: datos diarios de temperatura máxima. La variable debe ser de tipo double
  • tmin: datos diarios de temperatura mínima. La variable debe ser de tipo double

A continuación se muestran las primeras cinco filas del dataset y los tipos de datos de cada una de las variables.

knitr::kable(climate[1:10,])
date station_id tmax tmin prcp
1961-01-01 87448 37.4 13.5 0.6
1961-01-02 87448 27.4 14.3 23.9
1961-01-03 87448 26.6 13.5 0.0
1961-01-04 87448 31.0 11.7 6.0
1961-01-05 87448 27.0 14.1 0.0
1961-01-06 87448 26.3 11.3 0.0
1961-01-07 87448 34.1 12.0 6.7
1961-01-08 87448 32.8 15.9 0.0
1961-01-09 87448 37.6 16.1 0.0
1961-01-10 87448 26.9 4.6 0.0

Los nombres de las variables son importantes y deben ser siempre los mismos ya que el modelo las reconocerá a partir de los mismos. Los nombres deben ser los siguientes:

  • date : corresponde a la fecha del día en formato Date. El formato de la fecha para facilitar el reconocimiento por parte de R es “YYYY-MM-DD”, es decir, el año expresado con cuatro dígitos y luego dos dígitos para el mes y dos para el día.
  • station_id: Identificador unívoco de cada una de las estaciones. Debe ser un número entero.
  • tmax: temperatura máxima diaria expresada en °C.
  • tmin: temperatura mínima diaria expresada en °C.
  • prcp: precipitación diaria expresada en mm.

El orden de las variables no es importante pero, como se mencionó, si se deben respetar los nombres de cada una. En el caso de faltantes, no se utiliza ningún valor específico para los NAs, sólo se debe dejar ese valor vacío. Este archivo tiene un formato largo, es decir, las estaciones se deben colocar una debajo de la otra.

La generación de series sobre estaciones requiere de dos funciones básicas local_fit y local_simulate. Independientemente de si se incluyen totales trimestrales en el modelo, siempre se utilizan esas dos funciones.

Series sintéticas estacionarias

Ajuste de los modelos

A continuación se mostrará como ajustar los cuatro modelos estadísticos para una sola estación meteorológica para generar series estacionarias donde cada realización es complemente independiente.

El ajuste del modelo local (en un punto) necesita de dos funciones: en una se define la configuración general del modelo y con la segunda se corre el modelo propiamente dicho.

Primero se crea un objeto con el control para el ajuste del simulador. Los argumentos son:

  • prcp_occurrence_threshold: umbral de precipitación para un día lluvioso. La OMM recomienda un umbral de 0.1 mm para considerar un día como lluvioso.
  • avbl_cores: cantidad de núcleos disponibles para la paralelización.
  • planar_crs_in_metric_coords: sistema de coordenadas planar.
control_fit <-gamwgen::local_fit_control(
  prcp_occurrence_threshold = 0.1, # Umbral para la definición de días húmedos
  avbl_cores = 6, # Cantidad de núcleos disponibles
  planar_crs_in_metric_coords = 22185) # Sistema de referencia espacial (en metros)

Luego se corre el ajuste para la estación meteorológica con la función local_calibrate. Los argumentos de la función son:

  • climate: datos meteorológicos observados para la estación
  • stations: metadatos de las estaciones meteorológicas
  • seasonal_covariates: datos agregados trimestrales. Si es NULL el ajuste será sin covariables y las series generadas serán estacionarias.
  • control: objeto de control
  • verbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza el ajuste de los cuatro modelos para cada una de 
# las estaciones. En este caso, por cuestiones de tiempo a cargar un objeto ya precalculado. 
# Si el usuario desea correrlo deberá ver la nota anterior.
gamgen_fit <- gamwgen::local_calibrate(climate = climate, # Registro histórico de variables meteorológicas
  stations = stations, # Estaciones meteorológicas 
  seasonal_covariates = NULL, # Totales trimestrales de precipitación
  control = control_fit, # Objeto de control
  verbose = FALSE) # Impresión de mensajes en la consola.

Para esta demostración, cargamos el objeto con el ajuste del modelo ya realizado.

# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('input_data/local/fit_local_unconditional.RData')) {
  fs::file_copy(system.file('/autorun/local', "fit_local_unconditional.RData",  package = "gamwgen"),
              new_path = 'input_data/local/fit_local_unconditional.RData')
}

# Cargamos el archivo recientemente creado
load('input_data/local/fit_local_unconditional.RData')

# Clase del objeto con el ajuste del generador
class(gamgen_fit)
## [1] "gamwgen"
# Contenido del modelo 
names(gamgen_fit)
##  [1] "control"              "stations"             "climate"             
##  [4] "crs_used_to_fit"      "start_climatology"    "fitted_models"       
##  [7] "models_data"          "models_residuals"     "statistics_threshold"
## [10] "exec_times"

Dentro del objeto se guardan todo lo necesario para la simulación así como información accesoria.

  • control: copia de la configuración usada para calibrar el generador
  • stations: estaciones meteorológicas utilizadas para la calibración
  • climate: datos climáticos de cada uno de las estaciones
  • seasonal_covariates: series temporales de totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima.
  • crs_used_to_fit: sistema de referencia espacial usado para proyectar
  • start_climatology: climatología diaria de cada una de las variables de entrada.
  • fitted_models: modelos ajustados, uno para cada variable: temperaturas máxima y mínima y ocurrencia y montos de precipitación.
  • models_data: datos usados efectivamente usados para ajustar los modelos (sin NAs)
  • models_residuals: residuos de cada uno de los modelos. Es decir, la diferencia entre el valor ajustado por el modelo (clima local) y el valor observado en el día i
  • statistics_threshold: umbrales de amplitud térmica diaria por mes. Si la amplitud simulada está fuera de este rango, se repetirá la simulación para ese día a los fines de mantener la consistencia entre variables
  • exec_times: tiempo de ejecución de cada una de las etapas del ajuste

Cada uno de los GAMs ajustados se almacenan en el objeto gamgen_fit y pueden ser evaluados con la función summary().

summary(gamgen_fit$fitted_models$`87448`$tmax_fit)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## tmax ~ s(tmax_prev, tmin_prev, k = 50) + s(prcp_occ, bs = "re") + 
##     s(prcp_occ_prev, bs = "re") + s(doy, bs = "cc", k = 30)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.61343    0.03214   796.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df      F p-value    
## s(tmax_prev,tmin_prev) 30.2469  38.71  247.1  <2e-16 ***
## s(prcp_occ)             0.9989   1.00 1005.8  <2e-16 ***
## s(prcp_occ_prev)        0.9995   1.00 2300.2  <2e-16 ***
## s(doy)                 12.8918  28.00  158.9  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.71   Deviance explained = 71.1%
## fREML =  58955  Scale est. = 14.138    n = 21466

La función summary permite analizar los resultados del ajuste de cada uno de los modelos. Para el caso del modelo de temperatura máxima podemos ver la fórmula del GAM en la parte superior bajo el apartado Formula y la significancia de cada uno de los términos del modelo en la tabla inmediatamente inferior. Se puede observar que todos los términos son altamente significativos. También se incluyen como pruebas de bondad del ajuste el porcentaje de la varianza explicada por el modelo y el valor de R-ajustado.

El paquete gratia (Simpson, Gavin (2018)) es muy útil para la visualización de los ajustes de manera gráfica. Esta función toma los diagnósticos por defecto de la función gam.check del paquete mgcv (Wood, Simon (2001)) y los convierte en gráficos de la librería ggplot2 (Wickham, Hadley(2011)) que son visualmente muy atractivos. La figura está dividida en cuatro paneles, cada uno mostrando un diagnóstico diferente. En el panel superior izquierdo se muestra un Q-Q Plot de los residuos. Si el modelo ha ajustado satisfactoriamente los datos, los puntos deberían estar sobre la recta 1:1 demarcada en rojo. El panel superior derecho muestra los residuos vs el predictor lineal. El objetivo de este diagnóstico es que los residuos se distribuyan al azar y que no tengan un patrón claro. La presencia de patrones en los residuos indicaría que el modelo no ha explicado algun componente importante de la variabilidad de los datos. En el panel inferior izquierdo se muestra un histograma de los residuos siendo deseable que lo mismos tengan una distribución proxima a la Normal. El último gráfico en el panel inferior derecho muestra una gráfico de dispersión entre los valores ajustados y observados.

gratia::appraise(gamgen_fit$fitted_models$`87448`$tmax_fit) +
  ggplot2::theme_bw()

Estos diagnóstico exploratorios puede aplicarse a cada uno de los modelos ajustados por la función local_calibrate.

Con la creación del objeto gamgen_fit finaliza el proceso de ajuste y ya se pueden generar series sintéticas para la o las estaciones usadas para calibrar el generador.

Generación de series

La generación de series sigue la misma estructura anterior, una función para configurar la generación y otra que realiza la generación propiamente dicha.

Los argumentos de la función de control son:

  • nsim: cantidad de simulaciones a realizar. Se debe ingresar un valor entero mayor o igual a 1.
  • seed: semilla. Se debe ingresar cualquier numero entero. No es necesario recordarlo porque se guarda junto a los resultados.
  • avbl_cores: cantidad de núcleos disponibles para la paralelización.
  • use_spatially_correlated_noise: utilizar la generación estocástica espacialmente correlacionada. Esta opción sólo es válida si en el ajuste y en la simulación se usaron más de cinco estaciones meteorológicas diferentes. Con un menor número no es posible calcular los variogramas necesarios para la generación de los campos aleatorios. Se debe introducir un boolean (TRUE or FALSE).
  • use_temporary_files_to_save_ram: si se simulan muchas realizaciones o los recursos informáticos son escasos, esta opción permite guardar los resultados de cada una de las realizaciones en el disco liberando memoria RAM que quedará disponible para generar nuevas simulaciones. Al finalizar la generación todos los archivos se combinan en uno único. Se debe introducir un boolean (TRUE or FALSE).
  • use_temporary_files_to_save_ram: esta opción permite eliminar los archivos temporales creados para ahorrar RAM luego de terminar la generación de todas las simulaciones. Se debe introducir un boolean (TRUE or FALSE).
control_sim <- gamwgen::local_simulation_control(
  nsim = 10, # Cantidad de simulaciones a realizar
  seed = 1234, # Semilla para que los resultados sean reproducibles
  avbl_cores = 6, # Cantidad de núcleos disponibles a utilizar
  use_spatially_correlated_noise = FALSE, # Usar modelo de ruido espacialmente correlacionado
  use_temporary_files_to_save_ram = FALSE, # Guardar resultados intermedios para ahorrar RAM
  remove_temp_files_used_to_save_ram = TRUE) # Borrar los resultados intermedios creados anteriormente

Luego se procede a la simulación de datos meteorológicos. Los argumentos de la función de simulación son:

  • model: objeto con el resultado de la función local_calibrate()
  • simulation_locations: objeto tipo sf con la ubicación de las estaciones a simular. Las estaciones usadas deben haber sido incluidas en el proceso de ajuste. No es necesario que todas estén presentes, se pueden generar series solo sobre algunas de ellas.
  • start_date: fecha de comienzo de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de comienzo es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, ni tampoco posterior. Se debe introducir una fecha en formato date
  • end_date: fecha de fin de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de fin es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, tampoco puede ser posterior. Se debe introducir una fecha en formato date
  • control: objeto de control creado con la función control_sim().
  • output_folder: ruta al directorio donde se guardarán los resultados, tanto finales como intermedios.
  • output_filename: nombre del archivo de salida. Para facilitar la interoperabilidad, el archivo generado es un archivo de texto en formato separado por comas (.csv)
  • seasonal_covariates: datos agregados trimestrales. Si el ajuste se realizó con covariables, la generación también debe realizarse con ellas. Caso contrario se producirá un error. Se debe introducir un data frame con los valores agregados para las tres variables (precipitación y temperaturas máxima y mínima) pero no necesariamente deben ser los mismos a los utilizados en el ajuste. Si se desean simular tendencias de algún tipo, ya sea de un modelo de cambio climático o arbitrarias, se deben perturbar estas variables trimestrales e introducirlas aquí. Estas series si deben tener la misma longitud que el período a generar
  • verbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza la generación de series para cada una de las estaciones. 
# En este caso, por cuestiones de tiempo, vamos a cargar un objeto con los resultados de la simulación 
simulated_climate <- gamwgen::local_simulation(model = gamgen_fit, # Objeto con los resultados del ajuste
    simulation_locations = stations, # Estaciones para las cuales simular
    start_date = as.Date('2019-01-01'), # Fecha de comienzo de las simulaciones
    end_date = as.Date('2019-12-31'), # Fecha de fin de las simulaciones
    control = control_sim, # Objeto con la configuración
    output_folder = getwd(), # Directorio donde se guardarán los resultados
    output_filename = 'simulations.csv', # Nombre del archivo de salida
    seasonal_covariates = NULL, # Covariables estacionales
    verbose = FALSE) # Impresión de mensajes en la consola

Esta función produce dos tipos de resultados: una lista que permanece en el ambiente de R y los datos generados que son guardados como .csv en el directorio indicado precedentemente.

# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('output_data/simulated_local_unconditional.RData')) {
  fs::file_copy(system.file('/autorun/local', "simulated_local_unconditional.RData",  package = "gamwgen"),
              new_path = 'output_data/simulated_local_unconditional.RData')
}  
# Cargamos el archivo recientemente creado
load('output_data/simulated_local_unconditional.RData')

# Clase del objeto con el ajuste del generador
class(simulated_climate)
## [1] "list"            "gamwgen.climate"
# Contenido del modelo 
names(simulated_climate)
## [1] "nsim"                                       
## [2] "seed"                                       
## [3] "realizations_seeds"                         
## [4] "simulation_points"                          
## [5] "output_file_with_results"                   
## [6] "output_file_fomart"                         
## [7] "rdata_file_with_fitted_stations_and_climate"
## [8] "exec_times"

La lista contiene los siguientes objetos:

  • nsim: cantidad de realizaciones.
  • seed: semilla general para toda la generación. Corresponde a la que se incluye en la función de control.
  • realization_seeds: semillas para cada una de las realizaciones. Esto permite replicar los resultados.
  • simulation_points: puntos donde se generaron las series sintéticas.
  • output_fil_with_results: nombre del archivo con los resultados.
  • output_file_format: tipo de archivo de salida, en este caso .csv.
  • rdata_file_with_fitted_stations_and_climate: archivo .RData con los datos meteorológicos observados que fueron utilizados en el ajuste. También se incluyen los metadatos de cada uno de esos puntos.
  • exec_times: tiempo de ejecución del la generación.

Ahora veremos el formato del archivo de salida que contiene las series sintéticas.

Para desmenuzar la generación del tiempo local se pueden correr algunas de las funciones que forman parte de local_simulation. Por ejemplo, del objeto gamgen_fit se extraen los residuos y con la función gamwgen:::generate_residuals_statistics se calculan los parámetros.

En la tabla anterior se muestran los parámetros necesarios para generar el tiempo local de las temperaturas máxima y mínima.

gen_noise_params <- gamwgen:::generate_residuals_statistics(
            models_residuals = gamgen_fit$models_residuals)

#knitr::kable(gen_noise_params)
rmarkdown::paged_table(gen_noise_params)
  • station_id: número unívoco que identifica a cada estación meteorológica.
  • type: tipo de día lluvioso (Wet) o seco (Dry).
  • month: número de mes para los que se calculan los parámetros
  • sd.tmax_residuals: desvío estándar de los residuos del modelo de temperatura máxima.
  • sd.tmin_residuals: desvío estándar de los residuos del modelo de temperatura mínima.
  • mean.tmax_residuals: media de los residuos del modelo de temperatura máxima.
  • mean.tmin_residuals: media de los residuos del modelo de temperatura mínima * cov.residuals: covarianza de los residuos.
  • var.tmax_residuals: covarianza de los residuos del modelo de temperatura máxima.
  • var.tmin_residuals: covarianza de los residuos del modelo de temperatura mínima.

Losa parámetros para cada uno de los meses permiten generar valores de tiempo local para las dos temperaturas a partir de una distribución normal multivariada.

A continuación se muestra un ejemplo para el año 2019 sólo para los días lluviosos.

fechas <- data.frame(date = seq(as.Date('2019-01-01'), as.Date('2019-12-31'), 'days')) %>%
    dplyr::mutate(dia = lubridate::day(date),
           mes = lubridate::month(date))

# Ejemplo de generación de ruido para el mes de enero
tmax_dry <- purrr::map2_dfr(
    .x = fechas$dia,
    .y = fechas$mes,
    .f = function(dia, mes) {
        
        result_tmax_dry   <- control_sim$temperature_noise_generating_function(simulation_points =  stations %>%
                                                                               dplyr::filter(., station_id == '87448'),
                                                                           gen_noise_params = gen_noise_params,
                                                                           month_number = mes,
                                                                           selector = 'tmax_dry',
                                                                           seed = NULL)
        
        result_tmax_dry <- result_tmax_dry %>%
            dplyr::mutate(dia = dia, 
                          mes = mes)
        
    }
)  %>%
  sf::st_set_geometry(NULL) %>%
  dplyr::mutate(date = as.Date(paste0('2019-', mes, '-', dia))) %>%
  dplyr::select(-mes, -dia)

ggplot2::ggplot(data = tmax_dry  %>%
  tidyr::gather(residuo, valor, -date), 
  ggplot2::aes(x = date, y = valor, color = residuo)) +
  ggplot2::scale_y_continuous(limits = c(-15, 15), 
                              breaks = seq(-15, 15, 3), 
                              name = 'Tiempo local') +
  ggplot2::scale_x_date(name = 'Días') +
  ggplot2::scale_color_manual(values=c("DarkOrange", "Steelblue"),
                              labels = c("Temp. Máx.", "Temp. Min.")) +
  ggplot2::geom_line() +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom",
                 legend.title = ggplot2::element_blank())

La líonea naranja corresponde a la serie para temperaturas máximas y la azul para temperaturas mínimas.

# Se carga el set de datos simulados
simulated_climate <- readr::read_csv(here::here('output_data/local/simulated_local_unconditional.csv')) 
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   realization = col_double(),
##   station_id = col_double(),
##   point_id = col_double(),
##   longitude = col_double(),
##   latitude = col_double(),
##   date = col_date(format = ""),
##   tmax = col_double(),
##   tmin = col_double(),
##   prcp = col_double()
## )
# Primeras filas del objeto de salidas
knitr::kable(simulated_climate[1:10,])
realization station_id point_id longitude latitude date tmax tmin prcp
1 87448 1 5001636 6256577 2019-01-01 27.19266 5.464899 0
1 87448 1 5001636 6256577 2019-01-02 30.57296 10.266542 0
1 87448 1 5001636 6256577 2019-01-03 33.36446 15.180860 0
1 87448 1 5001636 6256577 2019-01-04 34.05627 15.473379 0
1 87448 1 5001636 6256577 2019-01-05 32.69240 15.270145 0
1 87448 1 5001636 6256577 2019-01-06 30.56683 14.797530 0
1 87448 1 5001636 6256577 2019-01-07 29.32913 12.624675 0
1 87448 1 5001636 6256577 2019-01-08 34.09988 13.816714 0
1 87448 1 5001636 6256577 2019-01-09 35.03660 19.596069 0
1 87448 1 5001636 6256577 2019-01-10 37.52834 19.323512 0

El resultado de la generación es un archivo .csv que contiene la siguiente información:

  • realization: número de realización. Es un valor entero entre 1 y la cantidad de realizaciones definida por el usuario.
  • station_id: número unívoco de identificación de la estación meteorológica o del punto arbitrario.
  • date: fechas de cada uno de los días de la simulación.
  • tmax: valores de temperatura máxima generada expresada en °C.
  • tmin: valores de temperatura mínima generada expresada en °C.
  • prcp: valores de precipitación diaria generada expresada en mm.

La siguiente Figura muestra un ejemplo de las series de temperaturas máximas y mínimas generadas.

# Grafico de temperatura máxima
tmax_plot <- ggplot2::ggplot() +
    ggplot2::geom_line(data = simulated_climate, 
                       ggplot2::aes(x = date, y = tmax, color = realization),
                       alpha = 0.5, color = 'DarkOrange') +
    ggplot2::geom_line(data = climate %>%
                           dplyr::filter(date > as.Date('2019-01-01')), 
                       ggplot2::aes(x = date, y = tmax)) +
    ggplot2::labs(x = 'Fecha', y = 'Temperatura máxima [°C]') +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "none")

# Grafico de temperatura mínima
tmin_plot <- ggplot2::ggplot() +
    ggplot2::geom_line(data = simulated_climate, 
                       ggplot2::aes(x = date, y = tmin, color = realization),
                       alpha = 0.5, color = 'Steelblue') +
    ggplot2::geom_line(data = climate %>%
                           dplyr::filter(date > as.Date('2019-01-01')), 
                       ggplot2::aes(x = date, y = tmin)) +
    ggplot2::labs(x = 'Fecha', y = 'Temperatura mínima [°C]') +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "none")

cowplot::ggdraw() +
    cowplot::draw_plot(tmax_plot, x = 0, y = 0.5, width = 1, height = .5) +
    cowplot::draw_plot(tmin_plot, x = 0, y = 0, width = 1, height = .5)

En el panel superior se muestra la temperatura máxima observada (negra) y las temperaturas sintéticas (naranja) para el año 2019. En el infeior se muestra la temperatura mínima diaria observada (negra) y cada una de las realizaciones (azul).

Series sintéticas pseudohistóricas

Ajuste de los modelos

A continuación se mostrará como ajustar los cuatro modelos estadísticos para una sola estación meteorológica para generar series pseudohistóricas donde cada realización copia las variaciones de baja frecuencia e la serie observada. El ajuste del modelo local (en un punto) necesita de dos funciones: en una se define la configuración general del modelo y con la segunda se corre el modelo propiamente dicho.

Primero se crea un objeto con el control para el ajuste del simulador. Los argumentos son:

  • prcp_occurrence_threshold: umbral de precipitación para un día lluvioso. La OMM recomienda un umbral de 0.1 mm para considerar un día como lluvioso.
  • avbl_cores: cantidad de núcleos disponibles para la paralelización.
  • planar_crs_in_metric_coords: sistema de coordenadas planar.
control_fit <-gamwgen::local_fit_control(
  prcp_occurrence_threshold = 0.1, # Umbral para la definición de días húmedos
  avbl_cores = 6, # Cantidad de núcleos disponibles
  planar_crs_in_metric_coords = 22185) # Sistema de referencia espacial (en metros)

Al tratarse de un modelo que ajusta condicionado por la variabilidad de baja frecuencia preexistente en los datos observados es necesario la agregación de las variables diarias en totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima. Esta operación puede realizarse con la función summarise_seasonal_climate incluida en el paquete. Esta función, además de agregar los datos, permite la imputación de faltantes. Se toleran una cierta cantidad que puede ser determinada por el usuario. El método de imputación utilizado es el imputePCA() de la librería missMDA .

# Agregación de valores diarios 
seasonal_covariates <- gamwgen::summarise_seasonal_climate(climate, umbral_faltantes = 0.2)

# Se muestran las primeras cinco filas
knitr::kable(seasonal_covariates[1:10,])
station_id year season seasonal_prcp seasonal_tmax seasonal_tmin
87448 1961 1 273.7 30.98764 14.5134831
87448 1961 2 236.0 24.47473 8.4670330
87448 1961 3 11.3 19.04565 1.3456522
87448 1961 4 234.5 24.66235 7.7611765
87448 1962 1 402.4 29.92333 14.2455556
87448 1962 2 187.2 24.30440 7.9868132
87448 1962 3 50.9 17.21957 -0.8978261
87448 1962 4 179.7 25.39560 7.8472527
87448 1963 1 306.8 29.09101 14.1404494
87448 1963 2 82.0 25.46196 8.5695652

Cabe mencionar que con esta función las valores se agregan por trimestre considerando la siguiente definición:

  • Verano: Diciembre, Enero y Febrero
  • Otoño: Marzo, Abril y Mayo
  • Invierno: Junio, Julio y Agosto
  • Primavera: Septiembre, Octubre y Noviembre

Los valores también se podrían agregar siguiendo otra definición de estaciones pero en ese caso, el usuario debería hacerlo por su cuenta. Algunas funciones útiles para hacerlo son las disponibles en el paquete lubridate como quarter() que permite definir el mes de comienzo de los trimestres. Para estas variables los nombres también son importantes por lo que deben respetarse los mostrados anteriormente.

La siguiente Figura muestra los totales trimestrales de precipitación para la estación meteorológica del ejemplo.

ggplot2::ggplot(data = seasonal_covariates %>%
                    dplyr::mutate(season = factor(season, 
                                                  labels = c('Verano', 'Otoño', 'Invierno', 'Primavera'))), 
                ggplot2::aes(x = year, y = seasonal_prcp,  group = 1)) +
    ggplot2::geom_line() +
    ggplot2::geom_smooth(se = FALSE, span = 0.5) +
    ggplot2::facet_wrap(~season, scales = 'free') +
    ggplot2::labs(x = 'Años', y = 'Precipitación acumulada [mm]') +
    ggplot2::theme_bw()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Cada uno de los paneles muestra la precipitación acumulada por trimestre junto a una regresión local (loess) en azul. Al incluir estos totales en el modelo, las series generadas tenderán a seguir la variabilidad observada, es decir, años secos generarán realización con valores nferiores a la media y viceversa.

Luego se corre el ajuste para la estación meteorológica con la función local_calibrate. Los argumentos de la función son:

  • climate: datos meteorológicos observados para la estación
  • stations: metadatos de las estaciones meteorológicas
  • seasonal_covariates: datos agregados trimestrales. Si es NULL el ajuste será sin covariables y las series generadas serán estacionarias.
  • control: objeto de control
  • verbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza el ajuste de los cuatro modelos para cada una de 
# las estaciones. En este caso, por cuestiones de tiempo a cargar un objeto ya precalculado. 
# Si el usuario desea correrlo deberá ver la nota anterior.
gamgen_fit <- gamwgen::local_calibrate(climate = climate, # Registro histórico de variables meteorológicas
  stations = stations, # Estaciones meteorológicas 
  seasonal_covariates = seasonal_covariates, # Totales trimestrales de precipitación
  control = control_fit, # Objeto de control
  verbose = FALSE) # Impresión de mensajes en la consola.

Para esta demostración, cargamos el objeto con el ajuste del modelo ya realizado.

# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('input_data/local/fit_local_conditional.RData')) {
  fs::file_copy(system.file('/autorun/local', "fit_local_conditional.RData",  package = "gamwgen"),
              new_path = 'input_data/local/fit_local_conditional.RData')
}

# Cargamos el archivo recientemente creado
load('input_data/local/fit_local_conditional.RData')

# Clase del objeto con el ajuste del generador
class(gamgen_fit)
## [1] "gamwgen"
# Contenido del modelo 
names(gamgen_fit)
##  [1] "control"              "stations"             "climate"             
##  [4] "seasonal_covariates"  "crs_used_to_fit"      "start_climatology"   
##  [7] "fitted_models"        "models_data"          "models_residuals"    
## [10] "statistics_threshold" "exec_times"

Dentro del objeto se guardan todo lo necesario para la simulación así como información accesoria.

  • control: copia de la configuración usada para calibrar el generador
  • stations: estaciones meteorológicas utilizadas para la calibración
  • climate: datos climáticos de cada uno de las estaciones
  • seasonal_covariates: series temporales de totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima.
  • crs_used_to_fit: sistema de referencia espacial usado para proyectar
  • start_climatology: climatología diaria de cada una de las variables de entrada.
  • fitted_models: modelos ajustados, uno para cada variable: temperaturas máxima y mínima y ocurrencia y montos de precipitación.
  • models_data: datos usados efectivamente usados para ajustar los modelos (sin NAs)
  • models_residuals: residuos de cada uno de los modelos. Es decir, la diferencia entre el valor ajustado por el modelo (clima local) y el valor observado en el día i
  • statistics_threshold: umbrales de amplitud térmica diaria por mes. Si la amplitud simulada está fuera de este rango, se repetirá la simulación para ese día a los fines de mantener la consistencia entre variables
  • exec_times: tiempo de ejecución de cada una de las etapas del ajuste

Cada uno de los GAMs ajustados se almacenan en el objeto gamgen_fit y pueden ser evaluados con la función summary().

summary(gamgen_fit$fitted_models$`87448`$tmax_fit)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## tmax ~ s(tmax_prev, tmin_prev, k = 50) + s(prcp_occ, bs = "re") + 
##     s(prcp_occ_prev, bs = "re") + s(doy, bs = "cc", k = 30) + 
##     s(SX1, SN1, k = 20) + s(SX2, SN2, k = 20) + s(SX3, SN3, k = 20) + 
##     s(SX4, SN4, k = 20)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.62116    0.03189   803.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df       F p-value    
## s(tmax_prev,tmin_prev) 29.9069 38.405  225.74  <2e-16 ***
## s(prcp_occ)             0.9989  1.000  930.30  <2e-16 ***
## s(prcp_occ_prev)        0.9995  1.000 2248.11  <2e-16 ***
## s(doy)                 12.6336 28.000   73.12  <2e-16 ***
## s(SX1,SN1)              3.0997  3.689   55.32  <2e-16 ***
## s(SX2,SN2)              2.0001  2.000   89.25  <2e-16 ***
## s(SX3,SN3)              4.4053  5.853   35.79  <2e-16 ***
## s(SX4,SN4)              2.0001  2.000   86.85  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.713   Deviance explained = 71.4%
## fREML =  58829  Scale est. = 13.964    n = 21466

La función summary permite analizar los resultados del ajuste de cada uno de los modelos. Para el caso del modelo de temperatura máxima podemos ver la fórmula del GAM en la parte superior bajo el apartado Formula y la significancia de cada uno de los términos del modelo en la tabla inmediatamente inferior. En este configuracion, al incluirle variables estacionales, se incoporan nuevos términos el modelo como SX1, SX2, SX3, SX4 y SN1, SN2, SN3, SN4, que corresponden a variables dummy de las medias trimestrales de temperatura máxima y mínima, respectivamente. Se puede observar que todos los términos son altamente significativos. También se incluyen como pruebas de bondad del ajuste el porcentaje de la varianza explicada por el modelo y el valor de R-ajustado.

El paquete gratia (Simpson, Gavin (2018)) es muy útil para la visualización de los ajustes de manera gráfica. Esta función toma los diagnósticos por defecto de la función gam.check del paquete mgcv (Wood, Simon (2001)) y los convierte en gráficos de la librería ggplot2 (Wickham, Hadley(2011)) que son visualmente muy atractivos. La figura está dividida en cuatro paneles, cada uno mostrando un diagnóstico diferente. En el panel superior izquierdo se muestra un Q-Q Plot de los residuos. Si el modelo ha ajustado satisfactoriamente los datos, los puntos deberían estar sobre la recta 1:1 demarcada en rojo. El panel superior derecho muestra los residuos vs el predictor lineal. El objetivo de este diagnóstico es que los residuos se distribuyan al azar y que no tengan un patrón claro. La presencia de patrones en los residuos indicaría que el modelo no ha explicado algun componente importante de la variabilidad de los datos. En el panel inferior izquierdo se muestra un histograma de los residuos siendo deseable que lo mismos tengan una distribución proxima a la Normal. El último gráfico en el panel inferior derecho muestra una gráfico de dispersión entre los valores ajustados y observados.

gratia::appraise(gamgen_fit$fitted_models$`87448`$tmax_fit) +
  ggplot2::theme_bw()

Estos diagnóstico exploratorios puede aplicarse a cada uno de los modelos ajustados por la función local_calibrate.

Con la creación del objeto gamgen_fit finaliza el proceso de ajuste y ya se pueden generar series sintéticas para la o las estaciones usadas para calibrar el generador.

Generación de series

La generación de series sigue la misma estructura anterior, una función para configurar la generación y otra que realiza la generación propiamente dicha.

Los argumentos de la función de control son:

  • nsim: cantidad de simulaciones a realizar. Se debe ingresar un valor entero mayor o igual a 1.
  • seed: semilla. Se debe ingresar cualquier numero entero. No es necesario recordarlo porque se guarda junto a los resultados.
  • avbl_cores: cantidad de núcleos disponibles para la paralelización.
  • use_spatially_correlated_noise: utilizar la generación estocástica espacialmente correlacionada. Esta opción sólo es válida si en el ajuste y en la simulación se usaron más de cinco estaciones meteorológicas diferentes. Con un menor número no es posible calcular los variogramas necesarios para la generación de los campos aleatorios. Se debe introducir un boolean (TRUE or FALSE).
  • use_temporary_files_to_save_ram: si se simulan muchas realizaciones o los recursos informáticos son escasos, esta opción permite guardar los resultados de cada una de las realizaciones en el disco liberando memoria RAM que quedará disponible para generar nuevas simulaciones. Al finalizar la generación todos los archivos se combinan en uno único. Se debe introducir un boolean (TRUE or FALSE).
  • use_temporary_files_to_save_ram: esta opción permite eliminar los archivos temporales creados para ahorrar RAM luego de terminar la generación de todas las simulaciones. Se debe introducir un boolean (TRUE or FALSE).
control_sim <- gamwgen::local_simulation_control(
  nsim = 10, # Cantidad de simulaciones a realizar
  seed = 1234, # Semilla para que los resultados sean reproducibles
  avbl_cores = 6, # Cantidad de núcleos disponibles a utilizar
  use_spatially_correlated_noise = FALSE, # Usar modelo de ruido espacialmente correlacionado
  use_temporary_files_to_save_ram = FALSE, # Guardar resultados intermedios para ahorrar RAM
  remove_temp_files_used_to_save_ram = TRUE) # Borrar los resultados intermedios creados anteriormente

Luego se procede a la simulación de datos meteorológicos. Los argumentos de la función de simulación son:

  • model: objeto con el resultado de la función local_calibrate()
  • simulation_locations: objeto tipo sf con la ubicación de las estaciones a simular. Las estaciones usadas deben haber sido incluidas en el proceso de ajuste. No es necesario que todas estén presentes, se pueden generar series solo sobre algunas de ellas.
  • start_date: fecha de comienzo de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de comienzo es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, ni tampoco posterior. Se debe introducir una fecha en formato date
  • end_date: fecha de fin de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de fin es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, tampoco puede ser posterior. Se debe introducir una fecha en formato date
  • control: objeto de control creado con la función control_sim().
  • output_folder: ruta al directorio donde se guardarán los resultados, tanto finales como intermedios.
  • output_filename: nombre del archivo de salida. Para facilitar la interoperabilidad, el archivo generado es un archivo de texto en formato separado por comas (.csv)
  • seasonal_covariates: datos agregados trimestrales. Si el ajuste se realizó con covariables, la generación también debe realizarse con ellas. Caso contrario se producirá un error. Se debe introducir un data frame con los valores agregados para las tres variables (precipitación y temperaturas máxima y mínima) pero no necesariamente deben ser los mismos a los utilizados en el ajuste. Si se desean simular tendencias de algún tipo, ya sea de un modelo de cambio climático o arbitrarias, se deben perturbar estas variables trimestrales e introducirlas aquí. Estas series si deben tener la misma longitud que el período a generar
  • verbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza la generación de series para cada una de las estaciones. 
# En este caso, por cuestiones de tiempo, vamos a cargar un objeto con los resultados de la simulación 
simulated_climate <- gamwgen::local_simulation(model = gamgen_fit, # Objeto con los resultados del ajuste
    simulation_locations = stations, # Estaciones para las cuales simular
    start_date = as.Date('2010-01-01'), # Fecha de comienzo de las simulaciones
    end_date = as.Date('2019-12-31'), # Fecha de fin de las simulaciones
    control = control_sim, # Objeto con la configuración
    output_folder = getwd(), # Directorio donde se guardarán los resultados
    output_filename = 'simulations.csv', # Nombre del archivo de salida
    seasonal_covariates = seasonal_covariates, # Covariables estacionales
    verbose = FALSE) # Impresión de mensajes en la consola

Esta función produce dos tipos de resultados: una lista que permanece en el ambiente de R y los datos generados que son guardados como .csv en el directorio indicado precedentemente.

# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('output_data/local/simulated_local_conditional.RData')) {
  fs::file_copy(system.file('/autorun/local', "simulated_local_conditional.RData",  package = "gamwgen"),
              new_path = 'output_data/local/simulated_local_conditional.RData')
}  
# Cargamos el archivo recientemente creado
load('output_data/local/simulated_local_conditional.RData')

# Clase del objeto con el ajuste del generador
class(simulated_climate)
## [1] "list"            "gamwgen.climate"
# Contenido del modelo 
names(simulated_climate)
## [1] "nsim"                                       
## [2] "seed"                                       
## [3] "realizations_seeds"                         
## [4] "simulation_points"                          
## [5] "output_file_with_results"                   
## [6] "output_file_fomart"                         
## [7] "rdata_file_with_fitted_stations_and_climate"
## [8] "exec_times"

La lista contiene los siguientes objetos:

  • nsim: cantidad de realizaciones.
  • seed: semilla general para toda la generación. Corresponde a la que se incluye en la función de control.
  • realization_seeds: semillas para cada una de las realizaciones. Esto permite replicar los resultados.
  • simulation_points: puntos donde se generaron las series sintéticas.
  • output_fil_with_results: nombre del archivo con los resultados.
  • output_file_format: tipo de archivo de salida, en este caso .csv.
  • rdata_file_with_fitted_stations_and_climate: archivo .RData con los datos meteorológicos observados que fueron utilizados en el ajuste. También se incluyen los metadatos de cada uno de esos puntos.
  • exec_times: tiempo de ejecución del la generación.

Ahora veremos el formato del archivo de salida que contiene las series sintéticas.

# Se carga el set de datos simulados
simulated_climate <- readr::read_csv(here::here('output_data/local/simulated_local_conditional.csv')) 
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   realization = col_double(),
##   station_id = col_double(),
##   point_id = col_double(),
##   longitude = col_double(),
##   latitude = col_double(),
##   date = col_date(format = ""),
##   tmax = col_double(),
##   tmin = col_double(),
##   prcp = col_double()
## )
# Primeras filas del objeto de salidas
knitr::kable(simulated_climate[1:10,])
realization station_id point_id longitude latitude date tmax tmin prcp
1 87448 1 5001636 6256577 2010-01-01 28.05965 17.69135 7.534622
1 87448 1 5001636 6256577 2010-01-02 27.32783 16.20019 0.000000
1 87448 1 5001636 6256577 2010-01-03 35.01025 13.35766 0.000000
1 87448 1 5001636 6256577 2010-01-04 34.18990 18.04710 0.000000
1 87448 1 5001636 6256577 2010-01-05 30.59661 20.62692 0.000000
1 87448 1 5001636 6256577 2010-01-06 37.15162 17.44571 0.000000
1 87448 1 5001636 6256577 2010-01-07 35.65166 18.36680 0.000000
1 87448 1 5001636 6256577 2010-01-08 28.28271 13.77686 0.000000
1 87448 1 5001636 6256577 2010-01-09 29.66045 13.73328 0.000000
1 87448 1 5001636 6256577 2010-01-10 28.07379 14.17866 11.043809

Al utilizar totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima, las series generadas capturan los variaciones de baja frecuencia que se observan en la serie histórica. A continuación se muestra una Figura que ilustra lo anterior para la temperatura máxima.

ggplot2::ggplot() +
    ggplot2::geom_boxplot(data = simulated_climate %>%
                              dplyr::mutate(month = lubridate::month(date),
                                            year = lubridate::year(date),
                                            date = as.Date(paste0(year, '-', month, '-', 15L))) , 
                          ggplot2::aes(x = date, y = tmax, group = date, fill = 'DarkOrange'), alpha = 0.1) +
    ggplot2::geom_line(data = climate %>%
                           dplyr::filter(date > as.Date('2010-01-01')) %>%
                           dplyr::mutate(month = lubridate::month(date),
                                         year = lubridate::year(date)) %>%
                           dplyr::group_by(year, month) %>%
                           dplyr::summarise(tmax = median(tmax, na.rm = TRUE)) %>%
                           dplyr::mutate(date = as.Date(paste0(year, '-', month, '-', 15L))),
                       ggplot2::aes(x = date, y = tmax, group = 1, color = 'DarkOrange')) +
    ggplot2::theme_bw() +
    ggplot2::labs(x = 'Fecha', y = 'Temperatura máxima [°C]') +
    ggplot2::scale_x_date(breaks = '6 months', labels = scales::date_format("%m-%y")) +
    ggplot2::scale_fill_discrete(name = "", labels = c("Simulado")) +
    ggplot2::scale_color_discrete(name = "", labels = c("Observado")) +
    ggplot2::theme(legend.position = 'bottom')
## `summarise()` regrouping output by 'year' (override with `.groups` argument)

Las cajas corresponde a las distintas realizaciones agregadas a escala mensual y la línea naranja corresponde a la temperatura media mensual calculada para los años 2010 a 2019. Se observa como las cajas suben y bajan al ritmo de la media observada y como capturan los pequeños cambios que ocurren en un año específico pero no en otros.

Series sintéticas correlacionadas espacialmente

Una tercera alternativa para la generación de datos sobre estaciones meteorológicas combina las antes mostradas pero generando el tiempo local con métodos que contemplan la autocorrelación espacial. Para utilizar esta alternativa se deben disponer de más de una sola estación porque de otro modo no se podrían calcular los parámetros de los variogramas necesarios para la generación del tiempo local. Mientras más estaciones haya mejor, pero si pueden generar campos espaciales confiables con alrededor de 10 puntos.

Crear archivos de entrada

Para este ejemplo se utilizan datos de varias estaciones pero el formato de los mismos es igual a los mostrados anteriormente.

Los archivos necesarios son:

  • stations.csv
  • climate.csv

Los datos meteorológicos se dividen en dos archivos separados: stations.csv y climate.csv. Los nombres de los mismos no deben ser necesariamente iguales a los usados aquí.

Los metadatos de las estaciones se alojan en el archivo stations.csv. Este archivo contiene la información de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:

  • station_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integer
  • latitude: latitud en grados decimales. La variable debe ser de tipo double
  • longitude: longitud en grados decimales. La variable debe ser de tipo double

La tabla puede tener más variables pero sólo se necesitan las anteriores.

A continuación se muestran la primera fila del dataset y los tipo de datos de cada una de las variables.

# Vista de los metadatos de la estación
knitr::kable(stations)
station_id nombre lat_dec lon_dec elev
86330 Artigas -30.4000 -56.5100 120.38
86350 Rivera -30.9000 -55.5400 241.94
86360 Salto -31.4300 -57.9800 41.00
86430 Paysandú -32.3500 -58.0400 61.12
86440 Melo -32.3700 -54.1900 100.36
86460 Paso de los Toros -32.8000 -56.5300 75.48
86490 Mercedes -33.2500 -58.0700 17.01
86560 Colonia -34.4500 -57.7700 18.00
86565 Rocha -34.4900 -54.3100 18.16
86580 Aeropuerto Carrasco -34.8300 -56.0100 32.88
90000001 Las Brujas -34.6719 -56.3400 32.00
90000002 La Estanzuela -34.3372 -57.6922 72.00
90000003 Tacuarembó -31.7089 -55.8267 115.00
90000004 Treinta y Tres -33.2750 -54.1722 22.00
90000005 Salto Grande -31.2728 -57.8908 47.00

El objeto stations debe ser convertido de tibble a sf. El sistema de referencia espacial debe ser planar. No es necesario un sistema de referencia espacial en particular, solamente las coordenadas deben estar expresadas en metros.

# Convertimos el objeto stations a sf y se convierte su proyección de WGS 1984 a
# UTM Zona 21 S
stations %<>% 
  sf::st_as_sf(coords = c("lon_dec", "lat_dec"), crs =  4326) %>%
  sf::st_transform(crs = 32721)

En el siguiente mapa muestra la distribución de las estaciones meteorológicas usadas en el ejemplo.

# Convertir el objeto con las estaciones a coordenadas geográficas 
stations_geo <- stations %>% 
    sf::st_transform(crs = 4326) %>%
    dplyr::mutate(lat = sf::st_coordinates(.)[,2],
                  lon = sf::st_coordinates(.)[,1])

# Creación del mapa con la distribución de las estaciones
leaflet::leaflet(data = stations_geo) %>%
    leaflet::addTiles() %>%
    leaflet::setView(lat = -33, lng = -56, zoom = 5) %>%
    leaflet::addCircleMarkers(lat = ~lat, lng = ~lon, radius = 3, 
                              fillColor = "#377eb8",
                              color = "#377eb8",
                              stroke = FALSE, fillOpacity = 1, opacity = 1,
                              popup = ~sprintf("<b>%s (%d)</b><br>Lat.: %.3f<br>Lon.: %.3f<br>Elev: %.0f m",
                                    nombre, station_id, lat, lon, elev)) 

La información climática se aloja en el archivo climate.csv. Este archivo contiene los datos de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:

  • date: fecha del dato. La variable debe ser de tipo date
  • station_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integer
  • prcp: datos diarios de precipitación La variable debe ser de tipo double
  • tmax: datos diarios de temperatura máxima. La variable debe ser de tipo double
  • tmin: datos diarios de temperatura mínima. La variable debe ser de tipo double

A continuación se muestran las primeras cinco filas del dataset y los tipos de datos de cada una de las variables.

knitr::kable(climate[1:10,])
date station_id tmax tmin prcp
1981-01-01 86330 34.6 21.7 0.0
1981-01-02 86330 33.8 21.2 19.5
1981-01-03 86330 28.3 19.4 0.0
1981-01-04 86330 30.3 17.4 0.0
1981-01-05 86330 33.4 17.4 21.0
1981-01-06 86330 32.8 20.4 8.0
1981-01-07 86330 30.6 17.9 0.0
1981-01-08 86330 30.2 19.1 0.0
1981-01-09 86330 31.3 20.2 0.0
1981-01-10 86330 30.2 20.6 0.0

En total se utilizan 15 estaciones meteorologicas, 10 provistas por INUMET (Instituyo Uruguayo de Mteorología) y 5 por INIA (Instituto Nacional de Investigación Agropecuaria).

unique(climate$station_id)
##  [1]    86330    86350    86360    86430    86440    86460    86490    86560
##  [9]    86565    86580 90000001 90000002 90000003 90000004 90000005

Ajuste de los modelos

El ajuste de los modelos es igual que para los dos ejemplos antes mostrados. Se utiliza la función control_fit para la congifuración del ajuste y local_calibrate para realizar el ajuste.

control_fit <-gamwgen::local_fit_control(
  prcp_occurrence_threshold = 0.1, # Umbral para la definición de días húmedos
  avbl_cores = 6, # Cantidad de núcleos disponibles
  planar_crs_in_metric_coords = 32721) # Sistema de referencia espacial (en metros)

Este ejemplo se realizará utilizando covriables trimestrales pero es válido para series estacionarias.

# Agregación de valores diarios 
seasonal_covariates <- gamwgen::summarise_seasonal_climate(climate, umbral_faltantes = 0.2)

# Se muestran las primeras cinco filas
knitr::kable(seasonal_covariates[1:10,])
station_id year season seasonal_prcp seasonal_tmax seasonal_tmin
86330 1981 1 458.9 30.32778 19.151111
86330 1981 2 370.4 26.44674 14.951087
86330 1981 3 211.2 19.76630 9.120652
86330 1981 4 272.0 24.73516 13.093407
86330 1982 1 465.9 30.33667 17.764444
86330 1982 2 229.0 26.74130 14.385870
86330 1982 3 365.4 19.48913 10.109783
86330 1982 4 651.7 24.68242 13.048352
86330 1983 1 609.0 31.39222 19.440000
86330 1983 2 471.2 23.80761 13.631522

En la siguiente Figura se muestra la variación de la precipitación acumulada estival para las 15 estaciones utilizadas.

ggplot2::ggplot(data = seasonal_covariates %>%
                    dplyr::filter(season == 1), 
                ggplot2::aes(x = year, y = seasonal_prcp)) +
    ggplot2::geom_line() +
    ggplot2::facet_wrap(~station_id, ncol = 3) +
    ggplot2::theme_bw() +
    ggplot2::labs(x = 'Años', y = 'Precipitación estival [mm]')

Se observa en la figura que si bien hay cierta coincidencia, no todos los picos y valles ocurren en el mismo momento. Utilizar un método que genere el tiempo local considerando estas variaciones espaciales es muy útil para representar feahacientemente los patrones espciales y temporales.

# Al correr la función se realiza el ajuste de los cuatro modelos para cada una de 
# las estaciones. En este caso, por cuestiones de tiempo a cargar un objeto ya precalculado. 
# Si el usuario desea correrlo deberá ver la nota anterior.
gamgen_fit <- gamwgen::local_calibrate(climate = climate, # Registro histórico de variables meteorológicas
  stations = stations, # Estaciones meteorológicas 
  seasonal_covariates = seasonal_covariates, # Totales trimestrales de precipitación
  control = control_fit, # Objeto de control
  verbose = FALSE) # Impresión de mensajes en la consola.

Para esta demostración, cargamos el objeto con el ajuste del modelo ya realizado.

# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('input_data/local/fit_local_correlated_weather.RData')) {
  fs::file_copy(system.file('/autorun/local', "fit_local_correlated_weather.RData",  package = "gamwgen"),
              new_path = 'input_data/local/fit_local_correlated_weather.Rdata')
}

# Cargamos el archivo recientemente creado
load('input_data/local/fit_local_correlated_weather.RData')

# Clase del objeto con el ajuste del generador
class(gamgen_fit)
## [1] "gamwgen"
# Contenido del modelo 
names(gamgen_fit)
##  [1] "control"              "stations"             "climate"             
##  [4] "seasonal_covariates"  "crs_used_to_fit"      "start_climatology"   
##  [7] "fitted_models"        "models_data"          "models_residuals"    
## [10] "statistics_threshold" "exec_times"

Dentro del objeto se guardan todo lo necesario para la simulación así como información accesoria.

  • control: copia de la configuración usada para calibrar el generador
  • stations: estaciones meteorológicas utilizadas para la calibración
  • climate: datos climáticos de cada uno de las estaciones
  • seasonal_covariates: series temporales de totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima.
  • crs_used_to_fit: sistema de referencia espacial usado para proyectar
  • start_climatology: climatología diaria de cada una de las variables de entrada.
  • fitted_models: modelos ajustados, uno para cada variable: temperaturas máxima y mínima y ocurrencia y montos de precipitación.
  • models_data: datos usados efectivamente usados para ajustar los modelos (sin NAs)
  • models_residuals: residuos de cada uno de los modelos. Es decir, la diferencia entre el valor ajustado por el modelo (clima local) y el valor observado en el día i
  • statistics_threshold: umbrales de amplitud térmica diaria por mes. Si la amplitud simulada está fuera de este rango, se repetirá la simulación para ese día a los fines de mantener la consistencia entre variables
  • exec_times: tiempo de ejecución de cada una de las etapas del ajuste Al inspeccionar el objeto con los resultados se puede ver que cada estación meteorológica tiene su propia sublista donde se almacenan los modelos para cada una de ellas.
names(gamgen_fit$fitted_models)
##  [1] "86330"    "86350"    "86360"    "86430"    "86440"    "86460"   
##  [7] "86490"    "86560"    "86565"    "86580"    "90000001" "90000002"
## [13] "90000003" "90000004" "90000005"

Cada uno de los GAMs ajustados se almacenan en el objeto gamgen_fit y pueden ser evaluados con la función summary().

summary(gamgen_fit$fitted_models$`86330`$tmax_fit)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## tmax ~ s(tmax_prev, tmin_prev, k = 50) + s(prcp_occ, bs = "re") + 
##     s(prcp_occ_prev, bs = "re") + s(doy, bs = "cc", k = 30) + 
##     s(SX1, SN1, k = 20) + s(SX2, SN2, k = 20) + s(SX3, SN3, k = 20) + 
##     s(SX4, SN4, k = 20)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26.18758    0.03397   770.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df      F p-value    
## s(tmax_prev,tmin_prev) 32.4286 40.704 235.57  <2e-16 ***
## s(prcp_occ)             0.9979  1.000 520.81  <2e-16 ***
## s(prcp_occ_prev)        0.9987  1.000 866.16  <2e-16 ***
## s(doy)                 10.3038 28.000  33.61  <2e-16 ***
## s(SX1,SN1)              3.6312  4.458  26.92  <2e-16 ***
## s(SX2,SN2)              2.0001  2.000  36.90  <2e-16 ***
## s(SX3,SN3)              4.3491  5.631  21.09  <2e-16 ***
## s(SX4,SN4)              2.0007  2.001  40.04  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.778   Deviance explained = 77.9%
## fREML =  35496  Scale est. = 8.8219    n = 14128

La función summary permite analizar los resultados del ajuste de cada uno de los modelos. Para el caso del modelo de temperatura máxima podemos ver la fórmula del GAM en la parte superior bajo el apartado Formula y la significancia de cada uno de los términos del modelo en la tabla inmediatamente inferior. En este configuracion, al incluirle variables estacionales, se incoporan nuevos términos el modelo como SX1, SX2, SX3, SX4 y SN1, SN2, SN3, SN4, que corresponden a variables dummy de las medias trimestrales de temperatura máxima y mínima, respectivamente. Se puede observar que todos los términos son altamente significativos. También se incluyen como pruebas de bondad del ajuste el porcentaje de la varianza explicada por el modelo y el valor de R-ajustado.

El paquete gratia (Simpson, Gavin (2018)) es muy útil para la visualización de los ajustes de manera gráfica. Esta función toma los diagnósticos por defecto de la función gam.check del paquete mgcv (Wood, Simon (2001)) y los convierte en gráficos de la librería ggplot2 (Wickham, Hadley(2011)) que son visualmente muy atractivos. La figura está dividida en cuatro paneles, cada uno mostrando un diagnóstico diferente. En el panel superior izquierdo se muestra un Q-Q Plot de los residuos. Si el modelo ha ajustado satisfactoriamente los datos, los puntos deberían estar sobre la recta 1:1 demarcada en rojo. El panel superior derecho muestra los residuos vs el predictor lineal. El objetivo de este diagnóstico es que los residuos se distribuyan al azar y que no tengan un patrón claro. La presencia de patrones en los residuos indicaría que el modelo no ha explicado algun componente importante de la variabilidad de los datos. En el panel inferior izquierdo se muestra un histograma de los residuos siendo deseable que lo mismos tengan una distribución proxima a la Normal. El último gráfico en el panel inferior derecho muestra una gráfico de dispersión entre los valores ajustados y observados.

gratia::appraise(gamgen_fit$fitted_models$`86330`$tmax_fit) +
  ggplot2::theme_bw()

Estos diagnóstico exploratorios puede aplicarse a cada uno de los modelos ajustados por la función local_calibrate.

Con la creación del objeto gamgen_fit finaliza el proceso de ajuste y ya se pueden generar series sintéticas para la o las estaciones usadas para calibrar el generador.

Generación de series

La generación de series sigue la misma estructura anterior, una función para configurar la generación y otra que realiza la generación propiamente dicha.

Los argumentos de la función de control son:

  • nsim: cantidad de simulaciones a realizar. Se debe ingresar un valor entero mayor o igual a 1.
  • seed: semilla. Se debe ingresar cualquier numero entero. No es necesario recordarlo porque se guarda junto a los resultados.
  • avbl_cores: cantidad de núcleos disponibles para la paralelización.
  • use_spatially_correlated_noise: utilizar la generación estocástica espacialmente correlacionada. Esta opción sólo es válida si en el ajuste y en la simulación se usaron más de cinco estaciones meteorológicas diferentes. Con un menor número no es posible calcular los variogramas necesarios para la generación de los campos aleatorios. Se debe introducir un boolean (TRUE or FALSE).
  • use_temporary_files_to_save_ram: si se simulan muchas realizaciones o los recursos informáticos son escasos, esta opción permite guardar los resultados de cada una de las realizaciones en el disco liberando memoria RAM que quedará disponible para generar nuevas simulaciones. Al finalizar la generación todos los archivos se combinan en uno único. Se debe introducir un boolean (TRUE or FALSE).
  • use_temporary_files_to_save_ram: esta opción permite eliminar los archivos temporales creados para ahorrar RAM luego de terminar la generación de todas las simulaciones. Se debe introducir un boolean (TRUE or FALSE).
control_sim <- gamwgen::local_simulation_control(
  nsim = 10, # Cantidad de simulaciones a realizar
  seed = 1234, # Semilla para que los resultados sean reproducibles
  avbl_cores = 6, # Cantidad de núcleos disponibles a utilizar
  use_spatially_correlated_noise = TRUE, # Usar modelo de ruido espacialmente correlacionado
  use_temporary_files_to_save_ram = TRUE, # Guardar resultados intermedios para ahorrar RAM
  remove_temp_files_used_to_save_ram = TRUE) # Borrar los resultados intermedios creados anteriormente

Luego se procede a la simulación de datos meteorológicos. Los argumentos de la función de simulación son:

  • model: objeto con el resultado de la función local_calibrate()
  • simulation_locations: objeto tipo sf con la ubicación de las estaciones a simular. Las estaciones usadas deben haber sido incluidas en el proceso de ajuste. No es necesario que todas estén presentes, se pueden generar series solo sobre algunas de ellas.
  • start_date: fecha de comienzo de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de comienzo es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, ni tampoco posterior. Se debe introducir una fecha en formato date
  • end_date: fecha de fin de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de fin es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, tampoco puede ser posterior. Se debe introducir una fecha en formato date
  • control: objeto de control creado con la función control_sim().
  • output_folder: ruta al directorio donde se guardarán los resultados, tanto finales como intermedios.
  • output_filename: nombre del archivo de salida. Para facilitar la interoperabilidad, el archivo generado es un archivo de texto en formato separado por comas (.csv)
  • seasonal_covariates: datos agregados trimestrales. Si el ajuste se realizó con covariables, la generación también debe realizarse con ellas. Caso contrario se producirá un error. Se debe introducir un data frame con los valores agregados para las tres variables (precipitación y temperaturas máxima y mínima) pero no necesariamente deben ser los mismos a los utilizados en el ajuste. Si se desean simular tendencias de algún tipo, ya sea de un modelo de cambio climático o arbitrarias, se deben perturbar estas variables trimestrales e introducirlas aquí. Estas series si deben tener la misma longitud que el período a generar
  • verbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza la generación de series para cada una de las estaciones. 
# En este caso, por cuestiones de tiempo, vamos a cargar un objeto con los resultados de la simulación 
simulated_climate <- gamwgen::local_simulation(model = gamgen_fit, # Objeto con los resultados del ajuste
    simulation_locations = stations, # Estaciones para las cuales simular
    start_date = as.Date('2019-01-01'), # Fecha de comienzo de las simulaciones
    end_date = as.Date('2019-12-31'), # Fecha de fin de las simulaciones
    control = control_sim, # Objeto con la configuración
    output_folder = getwd(), # Directorio donde se guardarán los resultados
    output_filename = 'simulations.csv', # Nombre del archivo de salida
    seasonal_covariates = seasonal_covariates, # Covariables estacionales
    verbose = FALSE) # Impresión de mensajes en la consola

Esta función produce dos tipos de resultados: una lista que permanece en el ambiente de R y los datos generados que son guardados como .csv en el directorio indicado precedentemente.

# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('output_data/local/simulatied_local_correlated_weather.RData')) {
  fs::file_copy(system.file('/autorun/local', "simulatied_local_correlated_weather.RData",  package = "gamwgen"),
              new_path = 'output_data/local/simulatied_local_correlated_weather.RData')
}  
# Cargamos el archivo recientemente creado
load('output_data/local/simulated_local_unconditional.RData')

# Clase del objeto con el ajuste del generador
class(simulated_climate)
## [1] "list"            "gamwgen.climate"
# Contenido del modelo 
names(simulated_climate)
## [1] "nsim"                                       
## [2] "seed"                                       
## [3] "realizations_seeds"                         
## [4] "simulation_points"                          
## [5] "output_file_with_results"                   
## [6] "output_file_fomart"                         
## [7] "rdata_file_with_fitted_stations_and_climate"
## [8] "exec_times"

La lista contiene los siguientes objetos:

  • nsim: cantidad de realizaciones.
  • seed: semilla general para toda la generación. Corresponde a la que se incluye en la función de control.
  • realization_seeds: semillas para cada una de las realizaciones. Esto permite replicar los resultados.
  • simulation_points: puntos donde se generaron las series sintéticas.
  • output_fil_with_results: nombre del archivo con los resultados.
  • output_file_format: tipo de archivo de salida, en este caso .csv.
  • rdata_file_with_fitted_stations_and_climate: archivo .RData con los datos meteorológicos observados que fueron utilizados en el ajuste. También se incluyen los metadatos de cada uno de esos puntos.
  • exec_times: tiempo de ejecución del la generación.

Ahora veremos el formato del archivo de salida que contiene las series sintéticas.

Para desmenuzar la generación del tiempo local se pueden correr algunas de las funciones que forman parte de local_simulation. Por ejemplo, del objeto gamgen_fit se extraen los residuos y con la función gamwgen:::generate_residuals_statistics se calculan los parámetros.

En la tabla anterior se muestran los parámetros necesarios para generar el tiempo local de las temperaturas máxima y mínima.

gen_noise_params <- gamwgen:::generate_month_params(
  residuals = gamgen_fit$models_residuals,
  observed_climate = gamgen_fit$models_data,
  stations = gamgen_fit$stations)

# Ejemplo de variograma de residuos de temperatura máxima para días secos
gen_noise_params[[1]]$variogram_parameters$tmax_dry  
## [1]   0.000000   6.978337 495.779099

Para cada una de las variables, precipitación y temperaturas máxima y mínima, se ajusta un variograma por máxima verosimilitud que permite representar la variabilidad espacial de los residuos de los modelos GAMs, que corresponden al tiempo local. En el variograma se muestran los tres parámetros nugget, psill y range. El rango de los datos es de 500 km, lo que es lógico si se considera que los puntos están distribuidamente en la República Oriental del Uruguay.

En la siguiente FIgura se muestra la variabilidad espacial del tiempo local de temperatura del día 1 de enero de 2019 para días secos. Cabe mencionar que se crean también los mismos datos para los días húmedos y en función de la ocurrencia de precipitación se selecciona uno u otro.

RandomFields::RFoptions(printlevel = 0, spConform = FALSE)

temp_local   <- control_sim$temperature_noise_generating_function(simulation_points =  stations %>%
                                                                           dplyr::mutate(latitude = sf::st_coordinates(.)[,2],
                                                                                         longitude = sf::st_coordinates(.)[,1]),
                                                                       gen_noise_params = gen_noise_params,
                                                                       month_number = 1L,
                                                                       selector = 'Dry',
                                                                       seed = NULL) %>%
     dplyr::mutate(date = as.Date(paste0('2019-', '01', '-', '01'))) %>%
     tidyr::gather(variable, valor, -c('date', 'geometry'))


# Mapa de Uruguay
uruguay <- raster::getData('GADM', country='URY', level=1) %>%
  sf::st_as_sf(.) %>%
  sf::st_transform(crs = 32721)

# Paleta de colores
colr <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(9, 'YlOrRd'))
# Quiebres de la escala de colores
quiebres <- c(-8.0, -6.0, -4.0,-2.0,  0,  2.0,  4.0, 6.0, 8.0)
# Etiquetas para los colores
etiquetas <- c('-8', '-6',  '-4', '-2',  '0',  '2',  '4',  '6', '8')
# Etiquetas de las facetas
labs <- c("Temp. Máx", "Temp. Min.")
names(labs) <- c("tmax_residuals", "tmin_residuals")

ggplot2::ggplot(data = temp_local) +
    ggplot2::geom_sf(ggplot2::aes(color = valor), size = 4) +
    ggplot2::scale_color_gradientn(name = "Tiempo local [°C]",
                         colours = colr(9),
                         breaks = quiebres,
                         labels = etiquetas,
                         limits = c(-8, 8)) +
    ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
    ggplot2::facet_wrap(~variable, labeller = ggplot2::labeller(variable = labs)) +
    ggplot2::theme_bw()

Se puede observar como el tiempo local tiene una estructura espacial que es representada a través del variograma y permite que los resultados para las distintas estaciones sean espacialmente consistentes. Esto quiere decir que estaciones meteorológicas cercanas tenderán a tener valores de tiempo local similares.

Generación de series sintéticas para una grilla regular

La generación de series sobre grillas regular es muy similar a lo mostrado anterioemente. Los fundamentos son los mismos, sólo los modelos estadísticos que modelan el clima local incoporan un nuevo término que permite incluir la dimensión espacial.

Crear archivos de entrada

Para este ejemplo se utilizan datos de varias estaciones pero el formato de los mismos es igual a los mostrados anteriormente.

Los archivos necesarios son:

  • stations.csv
  • climate.csv

Los datos meteorológicos se dividen en dos archivos separados: stations.csv y climate.csv. Los nombres de los mismos no deben ser necesariamente iguales a los usados aquí.

Los metadatos de las estaciones se alojan en el archivo stations.csv. Este archivo contiene la información de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:

  • station_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integer
  • latitude: latitud en grados decimales. La variable debe ser de tipo double
  • longitude: longitud en grados decimales. La variable debe ser de tipo double

La tabla puede tener más variables pero sólo se necesitan las anteriores.

A continuación se muestran la primera fila del dataset y los tipo de datos de cada una de las variables.

# Vista de los metadatos de la estación
knitr::kable(stations)
station_id nombre lat_dec lon_dec elev
86330 Artigas -30.4000 -56.5100 120.38
86350 Rivera -30.9000 -55.5400 241.94
86360 Salto -31.4300 -57.9800 41.00
86430 Paysandú -32.3500 -58.0400 61.12
86440 Melo -32.3700 -54.1900 100.36
86460 Paso de los Toros -32.8000 -56.5300 75.48
86490 Mercedes -33.2500 -58.0700 17.01
86560 Colonia -34.4500 -57.7700 18.00
86565 Rocha -34.4900 -54.3100 18.16
86580 Aeropuerto Carrasco -34.8300 -56.0100 32.88
90000001 Las Brujas -34.6719 -56.3400 32.00
90000002 La Estanzuela -34.3372 -57.6922 72.00
90000003 Tacuarembó -31.7089 -55.8267 115.00
90000004 Treinta y Tres -33.2750 -54.1722 22.00
90000005 Salto Grande -31.2728 -57.8908 47.00

El objeto stations debe ser convertido de tibble a sf. El sistema de referencia espacial debe ser planar. No es necesario un sistema de referencia espacial en particular, solamente las coordenadas deben estar expresadas en metros.

# Convertimos el objeto stations a sf y se convierte su proyección de WGS 1984 a
# UTM Zona 21 S
stations %<>% 
  sf::st_as_sf(coords = c("lon_dec", "lat_dec"), crs =  4326) %>%
  sf::st_transform(crs = 32721)

En el siguiente mapa muestra la distribución de las estaciones meteorológicas usadas en el ejemplo.

# Convertir el objeto con las estaciones a coordenadas geográficas 
stations_geo <- stations %>% 
    sf::st_transform(crs = 4326) %>%
    dplyr::mutate(lat = sf::st_coordinates(.)[,2],
                  lon = sf::st_coordinates(.)[,1])

# Creación del mapa con la distribución de las estaciones
leaflet::leaflet(data = stations_geo) %>%
    leaflet::addTiles() %>%
    leaflet::setView(lat = -33, lng = -56, zoom = 5) %>%
    leaflet::addCircleMarkers(lat = ~lat, lng = ~lon, radius = 3, 
                              fillColor = "#377eb8",
                              color = "#377eb8",
                              stroke = FALSE, fillOpacity = 1, opacity = 1,
                              popup = ~sprintf("<b>%s (%d)</b><br>Lat.: %.3f<br>Lon.: %.3f<br>Elev: %.0f m",
                                    nombre, station_id, lat, lon, elev)) 

La información climática se aloja en el archivo climate.csv. Este archivo contiene los datos de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:

  • date: fecha del dato. La variable debe ser de tipo date
  • station_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integer
  • prcp: datos diarios de precipitación La variable debe ser de tipo double
  • tmax: datos diarios de temperatura máxima. La variable debe ser de tipo double
  • tmin: datos diarios de temperatura mínima. La variable debe ser de tipo double

A continuación se muestran las primeras cinco filas del dataset y los tipos de datos de cada una de las variables.

knitr::kable(climate[1:10,])
date station_id tmax tmin prcp
1981-01-01 86330 34.6 21.7 0.0
1981-01-02 86330 33.8 21.2 19.5
1981-01-03 86330 28.3 19.4 0.0
1981-01-04 86330 30.3 17.4 0.0
1981-01-05 86330 33.4 17.4 21.0
1981-01-06 86330 32.8 20.4 8.0
1981-01-07 86330 30.6 17.9 0.0
1981-01-08 86330 30.2 19.1 0.0
1981-01-09 86330 31.3 20.2 0.0
1981-01-10 86330 30.2 20.6 0.0

En total se utilizan 15 estaciones meteorologicas, 10 provistas por INUMET (Instituyo Uruguayo de Meteorología) y 5 por INIA (Instituto Nacional de Investigación Agropecuaria).

unique(climate$station_id)
##  [1]    86330    86350    86360    86430    86440    86460    86490    86560
##  [9]    86565    86580 90000001 90000002 90000003 90000004 90000005

Ajuste de los modelos

El ajuste de los modelos es igual que para los dos ejemplos antes mostrados. Se utiliza la función control_fit para la congifuración del ajuste y spatial_calibrate para realizar el ajuste.

control_fit <-gamwgen::spatial_fit_control(
  prcp_occurrence_threshold = 0.1, # Umbral para la definición de días húmedos
  avbl_cores = 6, # Cantidad de núcleos disponibles
  planar_crs_in_metric_coords = 32721) # Sistema de referencia espacial (en metros)

Este ejemplo se realizará utilizando covriables trimestrales pero es válido para series estacionarias.

# Agregación de valores diarios 
seasonal_covariates <- gamwgen::summarise_seasonal_climate(climate, umbral_faltantes = 0.2)

# Se muestran las primeras cinco filas
knitr::kable(seasonal_covariates[1:10,])
station_id year season seasonal_prcp seasonal_tmax seasonal_tmin
86330 1981 1 458.9 30.32778 19.151111
86330 1981 2 370.4 26.44674 14.951087
86330 1981 3 211.2 19.76630 9.120652
86330 1981 4 272.0 24.73516 13.093407
86330 1982 1 465.9 30.33667 17.764444
86330 1982 2 229.0 26.74130 14.385870
86330 1982 3 365.4 19.48913 10.109783
86330 1982 4 651.7 24.68242 13.048352
86330 1983 1 609.0 31.39222 19.440000
86330 1983 2 471.2 23.80761 13.631522
# Al correr la función se realiza el ajuste de los cuatro modelos para cada una de 
# las estaciones. En este caso, por cuestiones de tiempo a cargar un objeto ya precalculado. 
# Si el usuario desea correrlo deberá ver la nota anterior.
gamgen_fit <- gamwgen::spatial_calibrate(climate = climate, # Registro histórico de variables meteorológicas
  stations = stations, # Estaciones meteorológicas 
  seasonal_covariates = seasonal_covariates, # Totales trimestrales de precipitación
  control = control_fit, # Objeto de control
  verbose = FALSE) # Impresión de mensajes en la consola.

Para esta demostración, cargamos el objeto con el ajuste del modelo ya realizado.

# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('input_data/spatial/fit_spatial_conditional.RData')) {
  fs::file_copy(system.file('/autorun/spatial', "fit_spatial_conditional.RData",  package = "gamwgen"),
              new_path = 'input_data/spatial/fit_spatial_conditional.Rdata')
}

# Cargamos el archivo recientemente creado
load('input_data/spatial/fit_spatial_conditional.RData')

# Clase del objeto con el ajuste del generador
class(gamgen_fit)
## [1] "gamwgen"
# Contenido del modelo 
names(gamgen_fit)
##  [1] "control"              "stations"             "climate"             
##  [4] "seasonal_covariates"  "crs_used_to_fit"      "start_climatology"   
##  [7] "fitted_models"        "models_data"          "models_residuals"    
## [10] "statistics_threshold" "exec_times"

Dentro del objeto se guardan todo lo necesario para la simulación así como información accesoria.

  • control: copia de la configuración usada para calibrar el generador
  • stations: estaciones meteorológicas utilizadas para la calibración
  • climate: datos climáticos de cada uno de las estaciones
  • seasonal_covariates: series temporales de totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima.
  • crs_used_to_fit: sistema de referencia espacial usado para proyectar
  • start_climatology: climatología diaria de cada una de las variables de entrada.
  • fitted_models: modelos ajustados, uno para cada variable: temperaturas máxima y mínima y ocurrencia y montos de precipitación.
  • models_data: datos usados efectivamente usados para ajustar los modelos (sin NAs)
  • models_residuals: residuos de cada uno de los modelos. Es decir, la diferencia entre el valor ajustado por el modelo (clima local) y el valor observado en el día i
  • statistics_threshold: umbrales de amplitud térmica diaria por mes. Si la amplitud simulada está fuera de este rango, se repetirá la simulación para ese día a los fines de mantener la consistencia entre variables
  • exec_times: tiempo de ejecución de cada una de las etapas del ajuste

A diferencia de la configuración que ajusta distintos modelos para cada estación meteorológica, en te caso se ajusta un sólo modelos para toda la región de interés. Al visualizar lo que está almacenado en la sublista fitted_models se ve que solo hay cuatro GAMs, dos para la precipitación y dos para las temperaturas máxima y mínima.

names(gamgen_fit$fitted_models)
## [1] "prcp_occ_fit" "prcp_amt_fit" "tmax_fit"     "tmin_fit"

Cada uno de los GAMs ajustados se almacenan en el objeto gamgen_fit y pueden ser evaluados con la función summary().

summary(gamgen_fit$fitted_models$tmax_fit)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## tmax ~ te(tmax_prev, tmin_prev, longitude, latitude, d = c(2, 
##     2), k = length(unique_stations)) + te(type_day, longitude, 
##     latitude, d = c(1, 2), bs = c("re", "tp"), k = length(unique_stations)) + 
##     te(type_day_prev, longitude, latitude, d = c(1, 2), bs = c("re", 
##         "tp"), k = length(unique_stations)) + te(doy, longitude, 
##     latitude, d = c(1, 2), bs = c("cc", "tp"), k = length(unique_stations)) + 
##     te(SX1, SN1, longitude, latitude, d = c(2, 2), k = length(unique_stations)) + 
##     te(SX2, SN2, longitude, latitude, d = c(2, 2), k = length(unique_stations)) + 
##     te(SX3, SN3, longitude, latitude, d = c(2, 2), k = length(unique_stations)) + 
##     te(SX4, SN4, longitude, latitude, d = c(2, 2), k = length(unique_stations))
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    25.91      29.52   0.878     0.38
## 
## Approximate significance of smooth terms:
##                                               edf  Ref.df       F p-value    
## te(tmax_prev,tmin_prev,longitude,latitude) 76.902  86.484 1267.40  <2e-16 ***
## te(type_day,longitude,latitude)            27.680  29.000  161.22  <2e-16 ***
## te(type_day_prev,longitude,latitude)       14.703  15.000  726.24  <2e-16 ***
## te(doy,longitude,latitude)                 47.907 195.000   77.44  <2e-16 ***
## te(SX1,SN1,longitude,latitude)             18.355  23.456   75.14  <2e-16 ***
## te(SX2,SN2,longitude,latitude)              6.005   6.009  184.32  <2e-16 ***
## te(SX3,SN3,longitude,latitude)             30.642  36.395   48.54  <2e-16 ***
## te(SX4,SN4,longitude,latitude)              6.002   6.004  193.99  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.779   Deviance explained = 77.9%
## fREML = 5.221e+05  Scale est. = 8.9846    n = 207269

La función summary permite analizar los resultados del ajuste de cada uno de los modelos. Para el caso del modelo de temperatura máxima podemos ver la fórmula del GAM en la parte superior bajo el apartado Formula y la significancia de cada uno de los términos del modelo en la tabla inmediatamente inferior. En este configuracion, al incluirle variables estacionales, se incoporan nuevos términos el modelo como SX1, SX2, SX3, SX4 y SN1, SN2, SN3, SN4, que corresponden a variables dummy de las medias trimestrales de temperatura máxima y mínima, respectivamente. Se puede observar que todos los términos son altamente significativos. También se incluyen como pruebas de bondad del ajuste el porcentaje de la varianza explicada por el modelo y el valor de R-ajustado.

El paquete gratia (Simpson, Gavin (2018)) es muy útil para la visualización de los ajustes de manera gráfica. Esta función toma los diagnósticos por defecto de la función gam.check del paquete mgcv (Wood, Simon (2001)) y los convierte en gráficos de la librería ggplot2 (Wickham, Hadley(2011)) que son visualmente muy atractivos. La figura está dividida en cuatro paneles, cada uno mostrando un diagnóstico diferente. En el panel superior izquierdo se muestra un Q-Q Plot de los residuos. Si el modelo ha ajustado satisfactoriamente los datos, los puntos deberían estar sobre la recta 1:1 demarcada en rojo. El panel superior derecho muestra los residuos vs el predictor lineal. El objetivo de este diagnóstico es que los residuos se distribuyan al azar y que no tengan un patrón claro. La presencia de patrones en los residuos indicaría que el modelo no ha explicado algun componente importante de la variabilidad de los datos. En el panel inferior izquierdo se muestra un histograma de los residuos siendo deseable que lo mismos tengan una distribución proxima a la Normal. El último gráfico en el panel inferior derecho muestra una gráfico de dispersión entre los valores ajustados y observados.

gratia::appraise(gamgen_fit$fitted_models$tmax_fit) +
  ggplot2::theme_bw()

Estos diagnóstico exploratorios puede aplicarse a cada uno de los modelos ajustados por la función spatial_calibrate.

Con la creación del objeto gamgen_fit finaliza el proceso de ajuste y ya se pueden generar series sintéticas para la o las estaciones usadas para calibrar el generador.

Generación de series

La generación de series sigue la misma estructura anterior, una función para configurar la generación y otra que realiza la generación propiamente dicha.

Los argumentos de la función de control son:

  • nsim: cantidad de simulaciones a realizar. Se debe ingresar un valor entero mayor o igual a 1.
  • seed: semilla. Se debe ingresar cualquier numero entero. No es necesario recordarlo porque se guarda junto a los resultados.
  • avbl_cores: cantidad de núcleos disponibles para la paralelización.
  • use_spatially_correlated_noise: utilizar la generación estocástica espacialmente correlacionada. Esta opción sólo es válida si en el ajuste y en la simulación se usaron más de cinco estaciones meteorológicas diferentes. Con un menor número no es posible calcular los variogramas necesarios para la generación de los campos aleatorios. Se debe introducir un boolean (TRUE or FALSE).
  • use_temporary_files_to_save_ram: si se simulan muchas realizaciones o los recursos informáticos son escasos, esta opción permite guardar los resultados de cada una de las realizaciones en el disco liberando memoria RAM que quedará disponible para generar nuevas simulaciones. Al finalizar la generación todos los archivos se combinan en uno único. Se debe introducir un boolean (TRUE or FALSE).
  • use_temporary_files_to_save_ram: esta opción permite eliminar los archivos temporales creados para ahorrar RAM luego de terminar la generación de todas las simulaciones. Se debe introducir un boolean (TRUE or FALSE).
control_sim <- gamwgen::spatial_simulation_control(
  nsim = 10, # Cantidad de simulaciones a realizar
  seed = 1234, # Semilla para que los resultados sean reproducibles
  avbl_cores = 6, # Cantidad de núcleos disponibles a utilizar
  use_spatially_correlated_noise = TRUE, # Usar modelo de ruido espacialmente correlacionado
  use_temporary_files_to_save_ram = TRUE, # Guardar resultados intermedios para ahorrar RAM
  remove_temp_files_used_to_save_ram = TRUE) # Borrar los resultados intermedios creados anteriormente
# Agregación de valores diarios 
seasonal_covariates <- gamwgen::summarise_seasonal_climate(climate, umbral_faltantes = 0.2)
# Shapefile de Uruguay
uruguay <- raster::getData('GADM', country='URY', level=1) %>%
    sf::st_as_sf(.) %>%
    sf::st_transform(crs = 32721)

# Creación de la grilla sobre la que se simulará
grilla_simulacion_centers <- uruguay %>%
  dplyr::select(geometry) %>%
    sf::st_transform(gamgen_fit$crs_used_to_fit) %>%
    sf::st_make_grid(cellsize = 25000, what = 'centers') %>%
    sf::st_sf(grid_id = 1:length(.), geometry = .) %>%
    sf::st_intersection(uruguay)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggplot2::ggplot() +
    ggplot2::geom_sf(data = grilla_simulacion_centers) +
    ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
    ggplot2::theme_bw()

Al utilizar covariables trimestrales es necesario que éstas se interpolen para cada uno de los puntos de la grilla. A continuación se muestran dos ejemplos, uno para precipitación acumulada para el trimestre estival de 2019 y la temperatura máxima media del mismo período.

simulation_dates <-
    tibble::tibble(date = seq.Date(from = as.Date('2019-01-01'),
                                   to = as.Date('2019-01-31'),
                                   by = "days")) %>%
    dplyr::mutate(year   = lubridate::year(date),
                  month  = lubridate::month(date),
                  season = lubridate::quarter(date, fiscal_start = 12))

grilla_simulacion_centers <- uruguay  %>%
    sf::st_transform(gamgen_fit$crs_used_to_fit) %>%
    sf::st_make_grid(cellsize = 25000, what = 'centers') %>%
    sf::st_sf(grid_id = 1:length(.), geometry = .) %>%
    sf::st_intersection(uruguay) %>%
    dplyr::select(geometry)  %>%
    dplyr::mutate(latitude = sf::st_coordinates(.)[,2],
                  longitude = sf::st_coordinates(.)[,1])
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
seasonal_covariates <-
    gamwgen:::get_covariates(model = gamgen_fit, 
                             simulation_points = grilla_simulacion_centers, 
                             seasonal_covariates, 
                             simulation_dates, 
                             control = control_sim)
## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(d$data): CRS object has comment, which is lost in output
## Warning in proj4string(newdata): CRS object has comment, which is lost in output

## Warning in proj4string(newdata): CRS object has comment, which is lost in output
## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(d$data): CRS object has comment, which is lost in output
## Warning in proj4string(newdata): CRS object has comment, which is lost in output

## Warning in proj4string(newdata): CRS object has comment, which is lost in output
## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output
## Warning in proj4string(d$data): CRS object has comment, which is lost in output
## Warning in proj4string(newdata): CRS object has comment, which is lost in output

## Warning in proj4string(newdata): CRS object has comment, which is lost in output
aux <- gamwgen:::sf2raster(seasonal_covariates %>%
                        dplyr::filter(year == 2019, season == 1),
                    variable = 'seasonal_prcp') %>%
    raster::as.data.frame(., xy = TRUE) %>%
    dplyr::rename(., 'longitude' = 'x', 'latitude' = 'y') %>%
    dplyr::mutate(., x = longitude, y = latitude) %>%
    sf::st_as_sf(., coords = c( 'x', 'y')) %>%
    sf::st_set_crs(., 32721)  %>%
    sf::st_intersection(uruguay) %>%
    dplyr::select(latitude, longitude, z, geometry)  
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
colores <- colorRampPalette(c("#cccccc", "#7bccc4", "#2b8cbe", "#0868ac", 
                              "#084081", "#807dba", "#6a51a3", "#54278f", "#3f007d", "#e7298a", 
                              "#ce1256", "#67001f"))

quiebres <- c(0.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0)

etiquetas <- c('0', '5', '10', '20', '30', '40', '50', '60', '70', '80', '90', '100', '110', '120', '130', '140', '150', '160')

ggplot(data = aux, 
             ggplot2::aes(x = longitude, y = latitude, fill = z)) +
    geom_tile(colour="black") +
    scale_fill_gradientn(name = "Precipitation [mm]", 
                         colours = colores(17), 
                         breaks = quiebres,
                         labels = etiquetas) +
    ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
    ggplot2::theme_bw() +
    ggplot2::theme(panel.background = element_rect(fill = "ivory"),
          axis.text = element_text(size = 11, colour = 1),
          legend.key.height = unit(2.5, "cm"))

aux <- gamwgen:::sf2raster(seasonal_covariates %>%
                               dplyr::filter(year == 2019, season == 1),
                           variable = 'seasonal_tmax') %>%
    raster::as.data.frame(., xy = TRUE) %>%
    dplyr::rename(., 'longitude' = 'x', 'latitude' = 'y') %>%
    dplyr::mutate(., x = longitude, y = latitude) %>%
    sf::st_as_sf(., coords = c( 'x', 'y')) %>%
    sf::st_set_crs(., 32721)  %>%
    sf::st_intersection(uruguay) %>%
    dplyr::select(latitude, longitude, z, geometry)  
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
colr <- colorRampPalette(RColorBrewer::brewer.pal(9, 'YlOrRd'))

quiebres <- c(0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0)

etiquetas <- c('0', '5', '10', '15', '20', '25', '30', '35', '40', '45')

ggplot(data = aux, 
       ggplot2::aes(x = longitude, y = latitude, fill = z)) +
    geom_tile(colour="black") +
    scale_fill_gradientn(name = "Temperatura máxima [°C]", 
                         colours = colr(9), 
                         breaks = quiebres,
                         labels = etiquetas) +
    ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
    ggplot2::theme_bw() +
    ggplot2::  labs(title = "", x = "", y = "") +
    ggplot2::theme(panel.background = element_rect(fill = "ivory"),
                   axis.text = element_text(size = 11, colour = 1))

Luego se procede a la simulación de datos meteorológicos. Los argumentos de la función de simulación son:

  • model: objeto con el resultado de la función local_calibrate()
  • simulation_locations: objeto tipo sf con la ubicación de las estaciones a simular. Las estaciones usadas deben haber sido incluidas en el proceso de ajuste. No es necesario que todas estén presentes, se pueden generar series solo sobre algunas de ellas.
  • start_date: fecha de comienzo de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de comienzo es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, ni tampoco posterior. Se debe introducir una fecha en formato date
  • end_date: fecha de fin de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de fin es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, tampoco puede ser posterior. Se debe introducir una fecha en formato date
  • control: objeto de control creado con la función control_sim().
  • output_folder: ruta al directorio donde se guardarán los resultados, tanto finales como intermedios.
  • output_filename: nombre del archivo de salida. Para facilitar la interoperabilidad, el archivo generado es un archivo de texto en formato separado por comas (.csv)
  • seasonal_covariates: datos agregados trimestrales. Si el ajuste se realizó con covariables, la generación también debe realizarse con ellas. Caso contrario se producirá un error. Se debe introducir un data frame con los valores agregados para las tres variables (precipitación y temperaturas máxima y mínima) pero no necesariamente deben ser los mismos a los utilizados en el ajuste. Si se desean simular tendencias de algún tipo, ya sea de un modelo de cambio climático o arbitrarias, se deben perturbar estas variables trimestrales e introducirlas aquí. Estas series si deben tener la misma longitud que el período a generar
  • verbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza la generación de series para cada una de las estaciones. 
# En este caso, por cuestiones de tiempo, vamos a cargar un objeto con los resultados de la simulación 
simulated_climate <- gamwgen::spatial_simulation(model = gamgen_fit, # Objeto con los resultados del ajuste
    simulation_locations = grilla_simulacion_centers, # Estaciones para las cuales simular
    start_date = as.Date('2019-01-01'), # Fecha de comienzo de las simulaciones
    end_date = as.Date('2019-01-31'), # Fecha de fin de las simulaciones
    control = control_sim, # Objeto con la configuración
    output_folder = getwd(), # Directorio donde se guardarán los resultados
    output_filename = 'simulations.csv', # Nombre del archivo de salida
    seasonal_covariates = seasonal_covariates, # Covariables estacionales
    verbose = FALSE) # Impresión de mensajes en la consola

Esta función produce dos tipos de resultados: una lista que permanece en el ambiente de R y los datos generados que son guardados como .csv en el directorio indicado precedentemente.

# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('output_data/spatial/simulated_spatial_conditional.RData')) {
  fs::file_copy(system.file('/autorun/spatial', "simulated_spatial_conditional.RData",  package = "gamwgen"),
              new_path = 'output_data/spatial/simulated_spatial_conditional.RData')
}  
# Cargamos el archivo recientemente creado
load('output_data/spatial/simulated_spatial_conditional.RData')

# Clase del objeto con el ajuste del generador
class(simulated_climate)
## [1] "list"            "gamwgen.climate"
# Contenido del modelo 
names(simulated_climate)
## [1] "nsim"                                       
## [2] "seed"                                       
## [3] "realizations_seeds"                         
## [4] "simulation_points"                          
## [5] "output_file_with_results"                   
## [6] "output_file_fomart"                         
## [7] "rdata_file_with_fitted_stations_and_climate"
## [8] "exec_times"

La lista contiene los siguientes objetos:

  • nsim: cantidad de realizaciones.
  • seed: semilla general para toda la generación. Corresponde a la que se incluye en la función de control.
  • realization_seeds: semillas para cada una de las realizaciones. Esto permite replicar los resultados.
  • simulation_points: puntos donde se generaron las series sintéticas.
  • output_fil_with_results: nombre del archivo con los resultados.
  • output_file_format: tipo de archivo de salida, en este caso .csv.
  • rdata_file_with_fitted_stations_and_climate: archivo .RData con los datos meteorológicos observados que fueron utilizados en el ajuste. También se incluyen los metadatos de cada uno de esos puntos.
  • exec_times: tiempo de ejecución del la generación.

Ahora veremos el formato del archivo de salida que contiene las series sintéticas.

Para desmenuzar la generación del tiempo local se pueden correr algunas de las funciones que forman parte de local_simulation. Por ejemplo, del objeto gamgen_fit se extraen los residuos y con la función gamwgen:::generate_residuals_statistics se calculan los parámetros.

En la tabla anterior se muestran los parámetros necesarios para generar el tiempo local de las temperaturas máxima y mínima.

gen_noise_params <- gamwgen:::generate_month_params(
  residuals = gamgen_fit$models_residuals,
  observed_climate = gamgen_fit$models_data,
  stations = gamgen_fit$stations)

# Ejemplo de variograma de residuos de temperatura máxima para días secos
gen_noise_params[[1]]$variogram_parameters$tmax_dry  
## [1]   0.000000   7.045242 493.298928

Explicar variograma.

Los parámetros para cada uno de los meses permiten generar valores de tiempo local para las dos temperaturas a partir de una distribución normal multivariada.

A continuación se muestra un ejemplo para el año 2019 sólo para los días lluviosos.

RandomFields::RFoptions(printlevel = 0, spConform = FALSE)

grilla_simulacion_centers <- uruguay  %>%
    sf::st_transform(gamgen_fit$crs_used_to_fit) %>%
    sf::st_make_grid(cellsize = 25000, what = 'centers') %>%
    sf::st_sf(grid_id = 1:length(.), geometry = .) %>%
    sf::st_intersection(uruguay) %>%
    dplyr::select(geometry)  %>%
    dplyr::mutate(latitude = sf::st_coordinates(.)[,2],
                  longitude = sf::st_coordinates(.)[,1])


# Mapa de Uruguay
uruguay <- raster::getData('GADM', country='URY', level=1) %>%
  sf::st_as_sf(.) %>%
  sf::st_transform(crs = 32721)

temp_local   <- control_sim$temperature_noise_generating_function(
    simulation_points =  grilla_simulacion_centers,
    gen_noise_params = gen_noise_params,
    month_number = 1L,
    selector = 'Dry',
    seed = 1234) %>%
    dplyr::mutate(date = as.Date(paste0('2019-', '01', '-', '01'))) %>%
    dplyr::filter(date == as.Date('2019-01-01')) %>%
    tidyr::gather(variable, valor, -c('geometry')) 

tmax_residuals <- gamwgen:::sf2raster(temp_local %>%
                                          dplyr::filter(variable == 'tmax_residuals'),
                                      'valor') %>%
    raster::as.data.frame(., xy = TRUE) %>%
    dplyr::rename(., 'longitude' = 'x', 'latitude' = 'y') %>%
    dplyr::mutate(., x = longitude, y = latitude) %>%
    sf::st_as_sf(., coords = c( 'x', 'y')) %>%
    sf::st_set_crs(., 32721)  %>%
    sf::st_intersection(uruguay) %>%
    dplyr::select(latitude, longitude, z, geometry)  


colr <- colorRampPalette(RColorBrewer::brewer.pal(9, 'YlOrRd'))

quiebres <- c(0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0)

etiquetas <- c('0', '5', '10', '15', '20', '25', '30', '35', '40', '45')

ggplot(data = tmax_residuals, 
       ggplot2::aes(x = longitude, y = latitude, fill = z)) +
    geom_tile(colour="black") +
    scale_fill_gradientn(name = "Temperatura máxima [°C]", 
                         colours = colr(9), 
                         breaks = quiebres,
                         labels = etiquetas) +
    ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
    ggplot2::theme_bw() +
    ggplot2::  labs(title = "", x = "", y = "") +
    ggplot2::theme(panel.background = element_rect(fill = "ivory"),
                   axis.text = element_text(size = 11, colour = 1))

La líonea naranja corresponde a la serie para temperaturas máximas y la azul para temperaturas mínimas.

# Se carga el set de datos simulados
simulated_climate <- readr::read_csv(here::here('output_data/spatial/simulated_spatial_conditional.csv')) 
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   realization = col_double(),
##   point_id = col_double(),
##   longitude = col_double(),
##   latitude = col_double(),
##   date = col_date(format = ""),
##   tmax = col_double(),
##   tmin = col_double(),
##   prcp = col_double()
## )
# Primeras filas del objeto de salidas
knitr::kable(simulated_climate[1:10,])
realization point_id longitude latitude date tmax tmin prcp
1 1 453759.4 6590418 2019-01-01 33.03589 19.72844 1.826571
1 2 553759.4 6590418 2019-01-01 31.77257 19.65962 27.883864
1 3 578759.4 6590418 2019-01-01 32.71455 20.47479 32.879767
1 4 428759.4 6615418 2019-01-01 30.64363 16.02136 7.032041
1 5 453759.4 6615418 2019-01-01 31.83107 18.12258 14.671730
1 6 478759.4 6615418 2019-01-01 33.32291 18.75998 14.781878
1 7 503759.4 6615418 2019-01-01 32.22893 18.78158 11.007127
1 8 528759.4 6615418 2019-01-01 31.95957 17.71497 24.732878
1 9 553759.4 6615418 2019-01-01 31.93247 17.65988 26.055372
1 10 453759.4 6640418 2019-01-01 30.67991 17.13227 16.806296

El resultado de la generación es un archivo .csv que contiene la siguiente información:

  • realization: número de realización. Es un valor entero entre 1 y la cantidad de realizaciones definida por el usuario.
  • station_id: número unívoco de identificación de la estación meteorológica o del punto arbitrario.
  • date: fechas de cada uno de los días de la simulación.
  • tmax: valores de temperatura máxima generada expresada en °C.
  • tmin: valores de temperatura mínima generada expresada en °C.
  • prcp: valores de precipitación diaria generada expresada en mm.

La siguiente Figura muestra un ejemplo de las series de temperaturas máximas y mínimas generadas.

uruguay <- raster::getData('GADM', country='URY', level=1) %>%
    sf::st_as_sf(.) %>%
    sf::st_transform(crs = 32721)

simulated_climate_geo <- simulated_climate %>%
    dplyr::filter(date == as.Date('2019-01-01'),
                  realization == 1) %>%
    sf::st_as_sf(coords = c('longitude', 'latitude'), crs = 32721) %>% 
    gamwgen:::sf2raster(., 'prcp') %>%
    raster::as.data.frame(., xy = TRUE) %>%
    dplyr::rename(., 'longitude' = 'x', 'latitude' = 'y') %>%
    dplyr::mutate(., x = longitude, y = latitude) %>%
    sf::st_as_sf(., coords = c( 'x', 'y')) %>%
    sf::st_set_crs(., 32721)  %>%
    sf::st_intersection(uruguay) %>%
    dplyr::select(latitude, longitude, z, geometry)  
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
colores <- colorRampPalette(c("#cccccc", "#7bccc4", "#2b8cbe", "#0868ac", 
                             "#084081", "#807dba", "#6a51a3", "#54278f", "#3f007d", "#e7298a", 
                             "#ce1256", "#67001f"))

quiebres <- c(0.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0)

etiquetas <- c('0', '5', '10', '20', '30', '40', '50', '60', '70', '80', '90', '100', '110', '120', '130', '140', '150', '160')


ggplot2::ggplot(data = simulated_climate_geo) +
    ggplot2::geom_tile(ggplot2::aes(x = longitude, y = latitude, fill = z)) +
    ggplot2::scale_fill_gradientn(name = "Precipitacópn [mm]",
                                   colours = colores(17),
                                   breaks = quiebres,
                                   labels = etiquetas,
                                   limits = c(0, 100)) +
    ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE)

#cowplot::ggdraw() +
#    cowplot::draw_plot(tmax_plot, x = 0, y = 0.5, width = 1, height = .5) +
#    cowplot::draw_plot(tmin_plot, x = 0, y = 0, width = 1, height = .5)

En el panel superior se muestra la temperatura máxima observada (negra) y las temperaturas sintéticas (naranja) para el año 2019. En el infeior se muestra la temperatura mínima diaria observada (negra) y cada una de las realizaciones (azul).

Agregar covariables interpoladas!!!!

Modelo no lineal
Modelos Lineales (Generalizados)

GLM: terminos lineales GAM: terminos no necesariamente lineales

Considerando \(n\) observaciones, \(x_{i}\), \(y_{i}\), donde \(y_{i}\) es una observación de una variable aleatoria, \(Y_{i}\), con una esperanza de \(\mu_{i}\)\(E(Y)\). Un modelo apropiado para expresar la relación entre \(x\) e \(y\) es el siguiente:

\[ y_i = \mu_i + \epsilon_i \] donde, \(\mu_i\) = \(x_{i}\beta\)

\(\beta\), parámetro desconocido, y \(\epsilon_i\) son variables aleatorias con media cero y cada una con variancia \(\sigma^2\). Por lo que el modelo establece que \(Y\) está dado por \(x\) multiplicada por una constante más un termino aleatorio. \(Y\) es un ejemplode variable respuesta mientras que \(x\) es un ejemplo de un predictor.

El modelo lineal se describe mediante la siguiente ecuación:

\[ y_i = \beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2 + \ldots + \epsilon_i \] \(y_i\) corresponde a la respuesta como una combinación lineal de covariables, \(x_{ji}\), con un offset.

La variable \(y_i\sim\) puede ser de cualquier distribución de la familia exponencial (Normal, Poisson, etc.)

El termino de error \(\epsilon_i\) está linealmente distribuido

Una vez obtenida la cantidad de simulaciones deseadas, se procede a la validación, que consiste en contrastar dichos valores con los registros observados para evaluar cuán bien éstas han capturado las propiedades estadísticas de los registros históricos. Por último, una vez que los resultados son validados, es posible simular series sintéticas con confianza tanto en localidades, con o sin datos históricos, como sobre una grilla regular que represente a una región.

El modelo es capaz de generar dos tipos de series sintéticas: no condicionadas y condicionadas. Las primeras, son series estacionarias, es decir, que no tienen tendencia temporal por lo que representan las secuencias medias de una localidad o región. Las series condicionadas, en cambio, incluyen en su proceso de ajuste variables auxiliares que reproducen cambios observados (ej., variabilidad climática natural) o proyectados en el clima. Estas variables auxiliares (también llamadas “covariables”) son totales trimestrales de lluvia acumulada y medias trimestrales de temperatura máxima y mínima. Su inclusión permite modelar cambios interdecadales o probabilidades estacionales, por lo que se pueden utilizar salidas de pronósticos estacionales, de modelos de cambio climático, o incluso valores arbitrarios definidos por el usuario para crear escenarios de interés (por ejemplo, un aumento sostenido de la precipitación). Si en cambio se utilizan covariables calculadas a partir de los datos, se obtienen series que hemos denominado pseudohistóricas.

Este es un ejemplo de aplicación del Generador Estocástico diario y multisitio basado en GAMs. Para ello se ajustará el modelo local del generador para una localidad de la Pampa Húmeda Argentina.

Is this linear?

set.seed(2) ## simulate some data...
dat <- gamSim(1, n = 400, dist = "normal", scale = 0.5, verbose=FALSE)
dat <- dat[,c("y", "x0", "x1", "x2", "x3")]
p <- ggplot(dat,aes(y=y,x=x1)) +
      geom_point() +
      theme_minimal()
print(p)

References

Kleiber, W., Katz, R.W., Rajagopalan, B., 2012. Daily spatiotemporal precipitation simulation using latent and transformed gaussian processes. Water Resources Research 48.

Kleiber, W., Katz, R.W., Rajagopalan, B., 2013. Daily minimum and maximum temperature simulation over complex terrain. Ann. Appl. Stat. 7, 588–612.

Simpson, G., 2018. Introducing gratia. From the Bottom of the Heap.

Verdin, A., 2016. Stochastic space-time modeling for agricultural decision support in the argentine pampas (Thesis).

Wickham, H., 2011. Ggplot2. Wiley Interdisciplinary Reviews: Computational Statistics 3, 180–185.

Wood, S.N., 2001. Mgcv: GAMs and generalized ridge regression for r. R news 1, 20–25.

Wood, S.N., 2017. Generalized additive models: An introduction with r. Chapman; Hall/CRC.